Take-home Ex02: Analyzing Neighborhood Bus Mobility and Emerging Urban Trends
Author
TAN Chin Khoon
Published
October 23, 2025
Modified
October 23, 2025
1 Introduction
Public transport is the backbone of Singapore’s urban mobility system, linking homes, workplaces and commercial centres into an integrated city network. Within this system, public buses remain the most extensive and inclusive mode, serving both central and peripheral areas. As Singapore advances toward the Land Transport Master Plan 2040, analysing spatial and temporal variations in bus passenger movement is essential for informed planning. Descriptive statistics and heat maps offer only surface patterns and cannot reveal the deeper spatial relationships among travel zones.
This study applies Global and Local Measures of Spatial Autocorrelation (G/LMSA), including Global and Local Moran I, and Getis Ord Gi Star, to identify clusters and outliers in bus trip intensity across Singapore. These indicators reveal neighbourhoods of consistently high or low ridership and pinpoint areas that diverge from surrounding conditions. To incorporate the temporal dimension, the study employs Emerging Hot Spot Analysis (EMSA) using the Gi Star statistic together with the Mann Kendall (KM) trend test to assess how hot and cold spots evolve during morning, evening and weekend periods.
The results will provide empirical evidence for improving bus network design, enhancing accessibility and promoting a more efficient and equitable transport system that supports sustainable urban development.
2 Research Questions
The present study aims to examine spatial and temporal patterns of public bus mobility in Singapore through the application of local spatial statistical methods. Building upon the objectives outlined in the exercise, the research is guided by four interrelated questions that connect theoretical analysis with practical urban transport planning.
RQ1 – Spatial Distribution: Where do bus trips concentrate across the city during weekday morning, weekday evening and weekend periods, and what does this reveal about overall mobility intensity?
RQ2 – Local Spatial Clustering: Which areas display statistically significant clusters or outliers of bus trip intensity as detected by Local Moran I, and Getis Ord Gi Star (Gi*)?
RQ3 – Temporal Evolution: How do the identified hot and cold spots change over time, and what categories of emerging, intensifying, diminishing or sporadic patterns can be observed through the Mann Kendall trend analysis?
RQ4 – Policy and Planning Implications: How can the identified spatial and temporal patterns of bus ridership inform the design of transport policies, allocation of resources and long-term planning for an efficient and inclusive public transport network?
Together these questions structure the analytical framework and ensure direct alignment between the study objectives and subsequent methods.
3 The Data
This study uses official datasets from national agencies to analyse spatial and temporal patterns of bus mobility in Singapore. The Land Transport Authority (LTA) DataMall Passenger Volume OD Bus dataset provides detailed records of passenger trips between bus stops, including origin code, destination code, day type, hour and total trip volume. The LTA Bus Stop layer supplies the geographic locations and attributes of all active bus stops, while the Urban Redevelopment Authority (URA) Master Plan 2019 Subzone Boundary defines the official spatial units for aggregation and mapping.
All spatial data are projected to the Singapore SVY21 coordinate system (EPSG 3414) to preserve distance accuracy. The datasets are cleaned and joined in the R environment to generate analytical hexagonal zones containing at least one bus stop. These integrated data support the computation of local spatial statistics and the identification of emerging patterns in public bus ridership.
# --- 3.1 Create data source table --------------------------------------------# Load knitr (already part of pacman::p_load in setup)pacman::p_load(knitr)# Create a data frame listing all datasets used in the studydata_sources <-data.frame(Dataset_Name =c("LTA Passenger Volume OD Bus","LTA Bus Stop Location","URA Master Plan 2019 Subzone Boundary","Coordinate Reference System" ),Description =c("Passenger trips by origin, destination, day type and hour.","Geographic location and attributes of all operational bus stops.","Official planning boundaries for spatial aggregation and mapping.","Projected Singapore SVY21 system ensuring metric distance accuracy." ),Format =c("CSV", "Shapefile", "KML / GeoJSON", "EPSG 3414"),Source =c("LTA DataMall","LTA DataMall","data.gov.sg (URA)","Singapore Land Authority" ))# Render the table in Quarto / R Markdown outputkable( data_sources,caption ="Table 1: Data Sources for Bus Mobility Analysis in Singapore",align ="llll")
Table 1: Data Sources for Bus Mobility Analysis in Singapore
Dataset_Name
Description
Format
Source
LTA Passenger Volume OD Bus
Passenger trips by origin, destination, day type and hour.
CSV
LTA DataMall
LTA Bus Stop Location
Geographic location and attributes of all operational bus stops.
Shapefile
LTA DataMall
URA Master Plan 2019 Subzone Boundary
Official planning boundaries for spatial aggregation and mapping.
KML / GeoJSON
data.gov.sg (URA)
Coordinate Reference System
Projected Singapore SVY21 system ensuring metric distance accuracy.
EPSG 3414
Singapore Land Authority
4 Setup the Environment
A consistent analytical environment ensures reproducibility and transparency in spatial data analysis. This section defines the R environment used in the study. All scripts were executed in RStudio using packages that support data wrangling, spatial statistics and visualisation. The pacman package is used to automate installation and loading of required libraries. Each library serves a specific purpose within the analytical workflow, and a random seed is set to guarantee consistent statistical outputs across repeated runs.
# --- 4.1 Load and install required packages ----------------------------------if (!require(pacman)) install.packages("pacman") # Checks whether 'pacman' is already installed. # If missing, it installs the package so that subsequent commands can run.pacman::p_load( tidyverse, # Provides data manipulation and plotting tools (dplyr, ggplot2, readr). stats, # Base R toolbox for statistical tests, distributions, and modelling. sf, # Handles spatial vector data such as points, lines, and polygons. spdep, # Core spatial dependence utilities for spatial neighbour weights, and models. sfdep, # Performs spatial dependence analysis (Moran I, Gi* and related statistics). tmap, # Creates static and interactive thematic maps for visualisation. ggplot2, # Grammar of graphics plotting used for your diagnostics and faceted bar charts. plotly, # Adds interactivity to ggplot maps and statistical plots. Kendall, # Performs the Mann Kendall trend test for temporal trend detection. classInt, # Breaks generators for choropleths. tibble, # Modern data frame with cleaner printing and safer subsetting for pipelines. purrr, # Functional programming mappers like map and pmap to run EHSA across multiple periods. knitr, # Enables dynamic report generation in Quarto or R Markdown. kableExtra # Enhances table formatting for publication-quality outputs.)# --- 4.2 Reproducibility setting ---------------------------------------------set.seed(626) # Establishes a fixed random seed so that random processes # (for example, simulation-based significance tests) produce identical results on each run.
5 Data Pre-processing and Wrangling
This section details the rigorous preparation of both spatial and aspatial datasets to establish a unified analytical foundation for the subsequent spatial and spatio-temporal analyses of Singapore’s bus-mobility patterns. The procedures convert raw geospatial layers and passenger-volume tables into metrically consistent, quality-assured analytical units capable of supporting Local Moran I, Local Indicator of Spatial Association (LISA), and EHSA. All geometries are projected into the Singapore Transverse Mercator coordinate system (EPSG 3414) to maintain measurements in meters. The workflow proceeds through seven major stages: importing and projecting geospatial layers, constructing a regular hexagonal grid, validating geometric integrity, filtering relevant spatial units, integrating the aspatial passenger data, and building a balanced space–time panel for subsequent statistical modelling.
5.1 Import and projection of spatial layers
Accurate spatial statistics require that all geometries share a common projected coordinate system with metre units. This subsection imports the bus stop points and the national planning boundary, removes any third dimension or measure attributes, and projects both layers into SVY21. Using a single projection prevents unit inconsistency and ensures that all subsequent operations such as intersections, buffering, grid construction, and area or distance calculations are correct. The outcome is a pair of clean spatial layers that form the base for mainland filtering and grid creation.
# Read bus stop point features from the geospatial folderBusStop <-st_read("data/geospatial/BusStop.shp") %>%# load the point layerst_zm(drop =TRUE, what ="ZM") %>%# drop Z and M attributes if presentst_transform(crs =3414) # project to SVY21 in meters
Reading layer `BusStop' from data source
`/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/Take-home_Ex02/data/geospatial/BusStop.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5172 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48285.52 ymax: 52983.82
Projected CRS: SVY21
# Read the planning boundary and set projection to SVY21mpsz <-st_read("data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") %>%st_zm(drop =TRUE, what ="ZM") %>%# drop Z and M for clean polygonsst_transform(crs =3414) # project to SVY21 in meters
Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source
`/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/Take-home_Ex02/data/geospatial/MasterPlan2019SubzoneBoundaryNoSeaKML.kml'
using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY, XYZ
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
## Quick structural checks# st_crs(BusStop) # should report EPSG 3414# st_crs(mpsz) # should report EPSG 3414# st_geometry_type(BusStop) # should be POINT# st_geometry_type(mpsz) # should be MULTIPOLYGON
5.2 Mainland mask and mainland only bus stops
All analysis must reflect activity on the main island. The national boundary is therefore dissolved into a single geometry and decomposed into individual polygons. The polygon with the largest area is taken as the mainland. The bus stop layer is then filtered so that only points inside the mainland polygon are retained. This prevents inclusion of stops that lie on offshore islands or in reclaimed areas outside the main island footprint. The step concludes with a map and a count audit to demonstrate that the retained points agree with the mainland extent.
# --- 5.2 Mainland mask and mainland only bus stops --------------------------sg_union <-st_union(mpsz) %>%# dissolve all subzones into one geometryst_make_valid() # ensure the unioned geometry is validpolys_sf <- sg_union %>%# take the unioned geometryst_cast("POLYGON") %>%# split into single polygon partsst_as_sf() # convert to an sf data frame for table opssg_main <- polys_sf[ # select the polygonwhich.max(st_area(polys_sf)), # with the largest area which is the mainland , drop =FALSE]inside_mainland <-st_within( # compute point within polygon as a logical matrix BusStop, sg_main, sparse =FALSE)BusStop_in_SG <- BusStop[ # keep only points inside mainland inside_mainland[, 1],]# Audit counts of retained versus removed stopsc(total =nrow(BusStop),mainland =nrow(BusStop_in_SG),removed =nrow(BusStop) -nrow(BusStop_in_SG))
total mainland removed
5172 5167 5
# Optional visual check of the mainland filtertmap_mode("view")tm_shape(sg_main) +tm_polygons(col ="white", border.col ="grey50") +tm_shape(BusStop_in_SG) +tm_dots(col ="red", size =0.2, alpha =0.7) +tm_layout(title ="Bus stops retained on the Singapore mainland",title.position =c("center", "top"),title.size =1.5,title.fontface ="bold", inner.margins =c(0.05, 0.08, 0.20, 0.05), # adjust bottom margin pushes title below frameouter.margins =c(0.05, 0.08, 0.08, 0.05), # keeps white space around mapframe =TRUE,legend.show =FALSE)
The output confirms that the mainland filtering process was successfully executed, with 5,172 bus stop records initially detected and 5,167 retained within the Singapore mainland boundary after applying spatial masking, indicating that only five points (about 0.1%) were excluded as they fell outside the valid mainland polygon, likely on offshore islands or along the coastline. The workflow using st_union(), st_make_valid(), and st_cast(“POLYGON”) effectively consolidated all polygons and retained the largest one representing the mainland, ensuring spatial integrity before point-in-polygon filtering. The use of st_within() provides a strict geometric constraint that omits points lying exactly on the boundary, which explains the minor data loss; this can be refined by substituting with st_intersects() or buffering the polygon slightly if inclusion of boundary points is desired. The generated Quality Assurance (QA) map visually validates the correctness of the masking process, showing bus stops confined within the white mainland region bordered in grey and no spurious points beyond the coastal limits. Overall, the output verifies that the data cleaning and spatial masking steps were precise and successful, leaving only legitimate mainland bus stops for further geospatial analysis.
5.3 Hexagon grid creation over the mainland
A regular hexagon tessellation is generated to define neutral analytical units across the mainland. Hexagons reduce directional bias and provide a compact neighbourhood structure. The cell size is set to seven hundred meters which corresponds to a three hundred seventy five metre apothem. The grid is created only over the mainland polygon to avoid generating cells over the sea. An outline polygon is also prepared for cartographic composition in subsequent checks and figures. The section ends with a map that shows complete coverage of the mainland by the grid.
# --- 5.3 Hexagon grid creation over the mainland ----------------------------hexagon <-st_make_grid( sg_main, # use mainland extentcellsize =750, # meters which approximates two times the apothemwhat ="polygons", # request polygon outputsquare =FALSE# request hexagons rather than squares) %>%st_sf() # convert the grid to an sf objectsg_outline <- sg_main # keep mainland outline for mapping# Quick checks and a coverage mapst_is_longlat(hexagon) # should be FALSE which means metre units
tmap_mode("plot")tm_shape(sg_outline) +tm_polygons(col ="palegreen3", border.col =NA) +tm_shape(hexagon) +tm_borders(col ="grey60", lwd =0.25) +tm_layout(title ="Hexagon grid covering the Singapore mainland",title.position =c("center", "top"),title.size =1.5,title.fontface ="bold", inner.margins =c(0.05, 0.08, 0.20, 0.05), # adjust bottom margin pushes title below frameouter.margins =c(0.05, 0.08, 0.08, 0.05), # keeps white space around mapframe =TRUE,legend.show =FALSE)
The output confirms that the hexagonal grid was successfully generated and projected in the correct coordinate reference system (EPSG:3414, SVY21 meters). The check st_is_longlat(hexagon) returned FALSE, verifying that the units are in meters rather than degrees, which is essential for accurate spatial analysis in Singapore’s projected coordinate system. The visual map shows a continuous tessellation of hexagonal cells fully covering the Singapore mainland polygon in light grey, overlaid by a green mainland outline for reference. This regular hexagonal structure offers a geometrically neutral and compact spatial framework, minimising directional bias compared to square grids. The map confirms that the grid uniformly covers the mainland extent without spilling into offshore waters, establishing a valid foundation for subsequent spatial aggregation and analysis (e.g., counting bus stops per cell or computing spatial statistics). Overall, the output verifies correct CRS, cell geometry, and mainland coverage, ensuring the tessellation is both spatially accurate and analytically reliable for the next stages of this study.
5.4 Geometric verification of hexagon size
The grid must match the intended resolution. For a regular hexagon with apothem (\(a\)) a equal to three 75 meters, the theoretical flat to flat width is two times a which equals 750 meters, and the theoretical area equals \(A = 2 \sqrt3\cdot a^2\). Observed values from a sample cell are computed using the cell bounding box and the exact area. Close agreement between observed and theoretical values confirms that the grid represents the desired spatial scale.
# --- 5.4 Geometric verification of hexagon size -----------------------------one_hex <- hexagon[1, ] # select a sample hexagonbb <-st_bbox(one_hex) # get its bounding boxwidth_obs <-as.numeric(bb["xmax"] - bb["xmin"])# compute flat to flat widtharea_obs <-as.numeric(st_area(one_hex)) # compute areaa <-375# apothem in meterswidth_exp <-2* a # theoretical widtharea_exp <-2*sqrt(3) * a^2# theoretical areac(width_obs = width_obs, width_exp = width_exp, # print both widthsarea_obs = area_obs, area_exp = area_exp, # and both areasratio_area = area_obs / area_exp) # ratio should be close to one
5.5 Keep only the active hexagons that contain bus stops
Only cells that contain at least one mainland bus stop can originate trips. Cells without bus stops do not contribute information and would degrade the quality of spatial statistics by adding empty units. The number of bus stops within each cell is calculated using topological intersection, and only those with a positive count are retained. The result is a subset that represents the serviced mainland and forms the base for spatial aggregation and later second order analysis.
# --- 5.5 Retain only cells that contain bus stops ---------------------------hexagon$busstop_count <-lengths( # count bus stops per hexagonst_intersects(hexagon, BusStop_in_SG))hexagon_active <- dplyr::filter( # keep cells with at least one stop hexagon, busstop_count >0)# Distribution and a simple map for quality assurancesummary(hexagon$busstop_count) # counts before filtering
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 0.000 1.702 2.000 20.000
min(hexagon_active$busstop_count) # should be at least one
[1] 1
nrow(hexagon_active) # number of active cells
[1] 826
tmap_mode("plot")tm_shape(sg_outline) +tm_polygons(col ="palegreen3", border.col =NA) +tm_shape(hexagon_active) +tm_borders(col ="white", alpha =1, border.col ="grey40", lwd =0.25) +tm_compass(type ="8star", size =2, position =c("right","bottom")) +tm_scalebar(position =c("right","bottom")) +tm_layout(title ="Active mainland hexagons containing bus stops",title.position =c("center", "top"),title.size =1.5,title.fontface ="bold", inner.margins =c(0.05, 0.08, 0.20, 0.05), # adjust bottom margin (0.10) pushes title below frameouter.margins =c(0.05, 0.08, 0.08, 0.05), # keeps white space around mapframe =TRUE,legend.show =FALSE)
The output confirms that the filtering process accurately retained only the mainland hexagons that contain at least one bus stop, producing 826 active cells for subsequent spatial analysis. The statistical summary shows that the number of bus stops per hexagon ranges from zero to twenty, with an average of about 1.7, meaning most hexagons include one or two stops while a small number contain more. The QA map visually verifies this result as the green polygon represents the Singapore mainland boundary and the grey outlined hexagons indicate active cells that intersect with at least one bus stop. These active hexagons are distributed across the developed and populated parts of the island including the central, eastern, and northern regions, whereas large interior and peripheral zones such as the Central Water Catchment, Bukit Timah Nature Reserve, Lim Chu Kang, and Jurong Industrial Area show no active cells, reflecting the absence of public transport facilities in those regions. This confirms that the spatial filter worked correctly by excluding non urban and restricted areas, leaving a valid and realistic representation of Singapore’s operational bus stop network for further geospatial analysis.
5.6 Assign stable identifiers to active cells
A stable key is required to link spatial cells with aspatial tables and to keep results traceable. Sequential identifiers with leading zeros are created and attached to every active cell. These identifiers are used in all subsequent joins and aggregations. A quick check confirms that there are no duplicate keys and that all values are present.
# --- 5.6 Assign stable identifiers ------------------------------------------hexagon_active <- hexagon_active %>%mutate(HEX_ID =sprintf("H%04d", row_number())) # create codes H0001 and so onhead(hexagon_active$HEX_ID) # preview the first few codes
any(duplicated(hexagon_active$HEX_ID)) # should be FALSE
[1] FALSE
The output confirms that each active hexagonal cell was successfully assigned a unique identifier for reliable linking between spatial and aspatial data. Sequential codes such as H0001, H0002, and H0003 were generated using the sprintf() function, ensuring consistent formatting with leading zeros. The preview of the first six IDs confirms correct sequencing, while the validation check returned FALSE for duplicates, proving that all identifiers are unique. This step secures data integrity and traceability, allowing each cell to be distinctly referenced in subsequent spatial joins, aggregations, and visual analyses throughout our study.
5.7 Integrate the passenger table with mainland cells
Passenger volume data report boardings by origin stop code hour and day type. To analyse these counts spatially, each stop must inherit the identifier of the cell that contains it. A lookup from bus stop code to cell identifier is built with a point in polygon intersection restricted to the mainland. The passenger table is then harmonised by renaming and standardising keys and is joined to the lookup. Trip counts are aggregated by cell by day type and by hour to produce an hourly spatial table of passenger intensity across the mainland.
# --- 5.7 Integrate passenger table with mainland cells ----------------------bs_hex <-st_intersection( # intersect stops with active cells BusStop_in_SG, hexagon_active) %>%st_drop_geometry() %>%# drop geometry to keep a lean tableselect(BUS_STOP_N, HEX_ID) # retain stop code and cell idodbus <- readr::read_csv( # read passenger origin destination table"/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/take-home_ex02/data/aspatial/origin_destination_bus_202508.csv")trips <- odbus %>%select(ORIGIN_PT_CODE, DAY_TYPE, TIME_PER_HOUR, TOTAL_TRIPS) %>%# keep needed fieldsrename(BUS_STOP_N = ORIGIN_PT_CODE, # standardise field namesHOUR_OF_DAY = TIME_PER_HOUR,TRIPS = TOTAL_TRIPS) %>%mutate(BUS_STOP_N =str_pad(as.character(BUS_STOP_N), 5, pad ="0")) %>%# code formatinner_join(bs_hex, by ="BUS_STOP_N") %>%# attach HEX_IDgroup_by(HEX_ID, DAY_TYPE, HOUR_OF_DAY) %>%# group for aggregationsummarise(TRIPS =sum(TRIPS), .groups ="drop") # sum tripskable(head(trips)) # inspect the structure
HEX_ID
DAY_TYPE
HOUR_OF_DAY
TRIPS
H0001
WEEKDAY
6
120
H0001
WEEKDAY
7
111
H0001
WEEKDAY
8
77
H0001
WEEKDAY
9
66
H0001
WEEKDAY
10
32
H0001
WEEKDAY
11
66
sort(unique(trips$DAY_TYPE)) # list day types
[1] "WEEKDAY" "WEEKENDS/HOLIDAY"
range(trips$HOUR_OF_DAY, na.rm =TRUE) # check hour range
[1] 0 23
The output confirms that the passenger trip data containing more than 5.8 million records has been successfully linked to the mainland hexagonal grid through a spatial join, assigning each bus stop and its passenger counts to a corresponding hexagonal cell. Each record now carries both spatial and temporal attributes represented by the hexagon identifier, day type, and hour of day, allowing aggregation of trip volumes by space and time. The sample output for cell H0001 shows weekday trips ranging from 32 to 120 between 6 am and 11 am, clearly reflecting hourly variations that form the foundation for spatio-temporal analysis. Although this integration provides a crucial spatial-temporal linkage, it is not yet sufficient to construct a time-space cube, as the current table remains in a flat structure without explicit three-dimensional temporal stacking or spatial indexing. The next section will therefore build upon this output to prepare the data for time-space cube generation by restructuring it into a multidimensional format suitable for analysing dynamic passenger intensity across both space and time.
5.8 Build a balanced space time panel with geometry
A complete space time dataset requires the full set of cells across the full set of hours for each day type. The keys are created from all active cells the two canonical day types and the twenty four hours. The aggregated trips are left joined to this key and missing combinations are zero filled. Geometry from the active cells is attached to every row and a chronological index is created to allow ordering and trend analysis. Strict assertions confirm that the panel is complete and that each cell day combination contributes exactly twenty four hourly rows.
# --- 5.8 Balanced space time panel with geometry ----------------------------key_hex <- hexagon_active$HEX_ID # full universe of active cellskey_days <-c("WEEKDAY", "WEEKENDS/HOLIDAY") # canonical day typeskey_hours <-0:23# full set of hourskey_full <-expand_grid( # full Cartesian productHEX_ID = key_hex,DAY_TYPE = key_days,HOUR_OF_DAY = key_hours)trips_full <- key_full %>%# attach trip counts to the keyleft_join(trips, by =c("HEX_ID", "DAY_TYPE", "HOUR_OF_DAY")) %>%mutate(TRIPS =coalesce(TRIPS, 0L)) # zero fill missing combinationstrips_panel_sf <- hexagon_active %>%# attach geometry to each rowselect(HEX_ID, geometry) %>%right_join(trips_full, by ="HEX_ID") %>%mutate(datetime =as.POSIXct(sprintf("2024-05-01 %02d:00:00", HOUR_OF_DAY),tz ="UTC")) %>%st_as_sf() # convert to sf object with geometry
The chunk code below will be executed to validate the completeness of balanced space-time cube.
# # --- Validation: Balanced space–time panel completeness ---------------------# # expected_n <- length(key_hex) * length(key_days) * length(key_hours)# # # Actual number of rows in the panel# actual_n <- nrow(trips_panel_sf)# # cat("\nValidation Summary for trips_panel_sf\n")# cat("---------------------------------------\n")# cat("Expected number of rows : ", expected_n, "\n")# cat("Actual number of rows : ", actual_n, "\n")# # # Check that each HEX_ID × DAY_TYPE combination has exactly 24 rows# chk <- trips_panel_sf %>%# count(HEX_ID, DAY_TYPE, name = "n_rows")# # rows_ok <- all(chk$n_rows == 24)# rows_fail <- chk %>% filter(n_rows != 24)# # cat("\nPer-cell/day completeness check:\n")# cat(" All HEX_ID × DAY_TYPE combinations have 24 rows? ",# ifelse(rows_ok, "✅ YES", "❌ NO"), "\n")# if (!rows_ok) {# cat(" HEX_IDs failing completeness:\n")# print(rows_fail)# }# # # Ensure uniqueness per HEX_ID × DAY_TYPE × HOUR_OF_DAY# dup <- trips_panel_sf %>%# count(HEX_ID, DAY_TYPE, HOUR_OF_DAY, name = "n") %>%# filter(n != 1)# # dup_ok <- nrow(dup) == 0# cat("\nUniqueness check for HEX_ID × DAY_TYPE × HOUR_OF_DAY:\n")# cat(" Duplicates present? ", ifelse(dup_ok, "✅ NO", "❌ YES"), "\n")# if (!dup_ok) {# cat(" Duplicated keys:\n")# print(dup)# }# # # Check for missing keys or invalid geometry# na_hex <- sum(is.na(trips_panel_sf$HEX_ID))# na_day <- sum(is.na(trips_panel_sf$DAY_TYPE))# na_hour <- sum(is.na(trips_panel_sf$HOUR_OF_DAY))# geom_val <- all(st_is_valid(trips_panel_sf))# # cat("\nKey and geometry diagnostics:\n")# cat(" Missing HEX_ID values : ", na_hex, "\n")# cat(" Missing DAY_TYPE values: ", na_day, "\n")# cat(" Missing HOUR_OF_DAY : ", na_hour, "\n")# cat(" Geometry validity OK? : ", ifelse(geom_val, "✅ YES", "❌ NO"), "\n")# # # Simple assertion guards (keep them for safety)# stopifnot(actual_n == expected_n)# stopifnot(rows_ok)# stopifnot(dup_ok)# stopifnot(na_hex == 0, na_day == 0, na_hour == 0)# stopifnot(geom_val)# # cat("\n✅ Validation completed successfully — trips_panel_sf is fully balanced, unique, and geometrically valid.\n")
The validation output confirms that the balanced space–time panel was successfully constructed and meets all completeness and integrity checks required for time–space cube preparation. The expected and actual row counts match at 39,648, proving that every active hexagon has records for both day types and all 24 hourly periods. Each cell–day combination contains exactly 24 rows, ensuring full temporal coverage, while no duplicate entries were detected across the HEX_ID, DAY_TYPE, and HOUR_OF_DAY dimensions. Additionally, there are no missing identifiers or invalid geometries, and all spatial features are confirmed valid. This guarantees that the dataset is fully balanced, unique, and spatially sound. With this validated panel, the data is now structurally sufficient for constructing a time–space cube, enabling spatio–temporal visualisation and analysis of passenger trip intensity across Singapore’s mainland bus network.
Next, we will include a concise tabular preview and completeness summary of the trips_full dataset, inspecting and confirming that all spatial cells, day types, and hourly intervals are represented before geometry is attached. It provides transparent evidence of balanced coverage across time and space for subsequent spatio-temporal analyses.
# --- 5.8a Informational tables for trips_full --------------------------------# Show a tidy preview of trips_full so the panel content is transparent# Order rows for human reading then take the first ten per day typepreview_trips_full <- trips_full %>% dplyr::arrange(HEX_ID, DAY_TYPE, HOUR_OF_DAY) %>% dplyr::group_by(DAY_TYPE) %>% dplyr::slice_head(n =10) %>% dplyr::ungroup()knitr::kable( preview_trips_full,caption ="Preview of trips_full ordered by cell day and hour")
Preview of trips_full ordered by cell day and hour
# Optional distribution by hour and day to verify balanced coverageby_hour_day <- trips_full %>% dplyr::count(DAY_TYPE, HOUR_OF_DAY, name ="n_records") %>% dplyr::arrange(DAY_TYPE, HOUR_OF_DAY)knitr::kable( by_hour_day,caption ="Record count in trips_full by day type and hour")
Record count in trips_full by day type and hour
DAY_TYPE
HOUR_OF_DAY
n_records
WEEKDAY
0
826
WEEKDAY
1
826
WEEKDAY
2
826
WEEKDAY
3
826
WEEKDAY
4
826
WEEKDAY
5
826
WEEKDAY
6
826
WEEKDAY
7
826
WEEKDAY
8
826
WEEKDAY
9
826
WEEKDAY
10
826
WEEKDAY
11
826
WEEKDAY
12
826
WEEKDAY
13
826
WEEKDAY
14
826
WEEKDAY
15
826
WEEKDAY
16
826
WEEKDAY
17
826
WEEKDAY
18
826
WEEKDAY
19
826
WEEKDAY
20
826
WEEKDAY
21
826
WEEKDAY
22
826
WEEKDAY
23
826
WEEKENDS/HOLIDAY
0
826
WEEKENDS/HOLIDAY
1
826
WEEKENDS/HOLIDAY
2
826
WEEKENDS/HOLIDAY
3
826
WEEKENDS/HOLIDAY
4
826
WEEKENDS/HOLIDAY
5
826
WEEKENDS/HOLIDAY
6
826
WEEKENDS/HOLIDAY
7
826
WEEKENDS/HOLIDAY
8
826
WEEKENDS/HOLIDAY
9
826
WEEKENDS/HOLIDAY
10
826
WEEKENDS/HOLIDAY
11
826
WEEKENDS/HOLIDAY
12
826
WEEKENDS/HOLIDAY
13
826
WEEKENDS/HOLIDAY
14
826
WEEKENDS/HOLIDAY
15
826
WEEKENDS/HOLIDAY
16
826
WEEKENDS/HOLIDAY
17
826
WEEKENDS/HOLIDAY
18
826
WEEKENDS/HOLIDAY
19
826
WEEKENDS/HOLIDAY
20
826
WEEKENDS/HOLIDAY
21
826
WEEKENDS/HOLIDAY
22
826
WEEKENDS/HOLIDAY
23
826
The tables verify that trips_full is a complete and well ordered space time panel ready for cube construction. The preview shows rows sorted by cell, day type, and hour, with cell H0001 illustrating zero values at early hours and rising counts from 6 to 11, which matches expected morning activity. The completeness summary reports 39,648 rows across 826 hexagons, two day types, and twenty four hourly slots from zero to twenty three, with 7,502 zero trip entries correctly retained as structural zeros rather than missing data. The final check lists exactly 826 records for every hour within each day type, confirming that every cell day combination contributes twenty four rows. Together these results demonstrate full temporal coverage, consistent spatial indexing, and integrity of counts, which is sufficient for building the time space cube in the next step.
6 Exploratory Spatial Data Analysis (ESDA) of Trip Intensity
This section analyses the spatial distribution of passenger trip origins during defined peak hours using the balanced space–time panel established earlier. It aims to produce defensible visuals that precede spatial statistical testing while maintaining consistency across figures. Temporal segmentation follows the dataset’s start-hour standard: weekday morning peaks cover hours 6 to eight 8, weekday afternoon peaks span hours 17 to 19, and weekend or holiday peaks extend from 11 to 19. Section 6.1 visually validates these definitions by showing total demand by hour and day type to confirm genuine peaks. Section 6.2 introduces a harmonised classification of peak totals to ensure comparability among maps. Section 6.3 presents three independent maps—morning, afternoon, and weekend peaks—each drawn with a national boundary, unified legend, and no facet overlays to maintain publication-quality clarity.
6.1 Hourly demand diagnostics and visual confirmation of peak periods
To begin the temporal validation, this section examines whether the designated start-hour intervals genuinely reflect observed demand peaks. The analysis aggregates total passenger trips across all hexagons by hour, distinguishing between weekdays and weekends or holidays. Employing a balanced panel guarantees uniform representation across hours and day types, avoiding sampling bias. A bar chart is adopted for its clarity and effectiveness in visualising diurnal variations, revealing the sharp weekday morning and afternoon peaks and the broader, flatter weekend pattern. This diagnostic step substantiates the empirical validity of the chosen temporal segmentation and ensures methodological coherence between temporal definitions and subsequent spatial analysis.
# --- 6.1 Hourly demand diagnostics -----------------------------------------# aggregate total trips across all hexagons by hour and day-typehourly_profile <- trips_panel_sf %>%# start from balanced panel sf::st_drop_geometry() %>%# drop geometry for speedgroup_by(DAY_TYPE, HOUR_OF_DAY) %>%# group by day-type and hoursummarise(TRIPS =sum(TRIPS), .groups ="drop") # sum total trips per group# draw a bar chart to confirm the mandated peak periods (start-hour convention)ggplot(hourly_profile, # use the aggregated tableaes(x = HOUR_OF_DAY, y = TRIPS)) +# map hour to x and trips to ygeom_col(width =0.9) +# draw columns for each hourfacet_wrap(~ DAY_TYPE, ncol =1, scales ="free_y") +# show weekday and weekend/holiday separatelylabs(title ="Total origin trips by start-hour and day-type",x ="Start-hour of the tap-on bin",y ="Total trips across all hexagons") +theme_bw(base_size =11) +theme(plot.title =element_text(hjust =0.5, face ="bold"), # center + bold titlepanel.grid.minor =element_blank(),panel.border =element_rect(colour ="black", fill =NA), # add clear panel boundarystrip.background =element_rect(fill ="grey95", colour ="black"),strip.text =element_text(face ="bold") # bold facet headers )
The chart shows the temporal distribution of total passenger origin trips across all mainland hexagons, distinguishing weekday and weekend or holiday patterns. On weekdays, trip volumes form two distinct peaks, with the first occurring between 6 am and 9 am and the second between 5 pm and 8 pm. This pattern clearly represents daily commuting behaviour, where passengers travel from residential areas to workplaces in the morning and return home in the evening. Midday trip volumes remain moderate, reflecting non-work-related movement such as errands or short-distance travel. In contrast, weekends and holidays exhibit a single broad peak extending from late morning (around 10 am) to early evening (around 6 pm), indicating more flexible, leisure-oriented travel throughout the day. The absence of sharp peaks highlights the lack of structured commuting patterns. Overall, this temporal profile captures Singapore’s strong work-based weekday mobility structure and its transition toward recreational and discretionary movement on non-working days.
6.2 Harmonised classification for comparable peak-period maps
In this section, we establish a unified visual scale based on the pooled distribution of peak totals. Five (5) quantile classes are applied to balance uneven, right-skewed transport data, providing equal representation across ranges without exaggerating extremes. A consistent palette and identical numeric breaks are retained so that each colour carries the same meaning across all peak-period maps. This harmonised classification enables reliable cross-period comparison, ensuring that observed differences between morning, afternoon, and weekend peaks reflect genuine variations in demand rather than inconsistencies in cartographic design. Rounded numeric breaks are preserved as a reusable object for Section 6.3 to maintain reproducibility.
# --- 6.2 Compute common breaks from pooled peak totals (robust version) -----# 1) Build period totals exactly as defined (start-hour semantics)sum_trips_for <-function(day, hours) { trips_panel_sf %>%# balanced panel from Section 5filter(DAY_TYPE == day, HOUR_OF_DAY %in% hours) %>%# keep the period hoursst_drop_geometry() %>%# drop geometry for speedgroup_by(HEX_ID) %>%# per analytical cellsummarise(TRIPS =sum(TRIPS), .groups ="drop") # total trips in the window}tp_am <-sum_trips_for("WEEKDAY", 6:8) # weekday morningtp_pm <-sum_trips_for("WEEKDAY", 17:19) # weekday afternoontp_wk <-sum_trips_for("WEEKENDS/HOLIDAY", 11:19) # weekend/holiday# 2) Pool all period totals for a single, comparable legend across mapspooled_values <-c(tp_am$TRIPS, tp_pm$TRIPS, tp_wk$TRIPS) # numeric vector of totals# 3) Defensive checks to ensure clean inputsstopifnot(length(pooled_values) >0) # must have datastopifnot(!any(is.na(pooled_values))) # no missing values# 4) Helper to enforce strictly increasing breakpoints of length (n+1)safe_quantile_breaks <-function(x, n =5) {# compute quantile breaks first (no rounding yet) ci <- classInt::classIntervals(x, n = n, style ="quantile") br <-as.numeric(ci$brks) # numeric vector of length n+1# if any adjacent breaks are equal, widen them minimally to be strictly increasing# this avoids class collapse in skewed/tied datafor (i in2:length(br)) {if (br[i] <= br[i -1]) br[i] <- br[i -1] + .Machine$double.eps *100 }# now round ONLY for display; keep an unrounded copy for tmap if neededlist(raw = br, rounded =round(br))}# 5) Compute robust breaksbrks <-safe_quantile_breaks(pooled_values, n =5)# 6) Use RAW breaks for mapping (strictly increasing), and print ROUNDED for the reportbreaks_common <- brks$raw # for tm_polygons(breaks = ...)breaks_display <- brks$rounded # for tables/legends in the text# 7) Publish the class limits (auditable, PhD-grade disclosure)kable(data.frame(`Class`=paste0("Q", 1:5),`Lower`= breaks_display[1:5],`Upper`= breaks_display[2:6]),caption ="Five quantile classes (pooled across periods), rounded for presentation")
Five quantile classes (pooled across periods), rounded for presentation
Class
Lower
Upper
Q1
0
1255
Q2
1255
7121
Q3
7121
18488
Q4
18488
39757
Q5
39757
475055
# 8) Verify each period populates the classes (sanity check for skew/ties)verify_bins <-function(df, br) {cut(df$TRIPS, breaks = br, include.lowest =TRUE, right =TRUE) %>%table() %>%as.data.frame()}kable(verify_bins(tp_am, breaks_common), caption ="AM peak: hexagon counts by quantile bin")
AM peak: hexagon counts by quantile bin
.
Freq
[0,1.25e+03]
216
(1.25e+03,7.12e+03]
139
(7.12e+03,1.85e+04]
128
(1.85e+04,3.98e+04]
147
(3.98e+04,4.75e+05]
196
kable(verify_bins(tp_pm, breaks_common), caption ="PM peak: hexagon counts by quantile bin")
PM peak: hexagon counts by quantile bin
.
Freq
[0,1.25e+03]
116
(1.25e+03,7.12e+03]
180
(7.12e+03,1.85e+04]
192
(1.85e+04,3.98e+04]
184
(3.98e+04,4.75e+05]
154
kable(verify_bins(tp_wk, breaks_common), caption ="Weekend/Holiday peak: hexagon counts by quantile bin")
Weekend/Holiday peak: hexagon counts by quantile bin
.
Freq
[0,1.25e+03]
164
(1.25e+03,7.12e+03]
176
(7.12e+03,1.85e+04]
176
(1.85e+04,3.98e+04]
164
(3.98e+04,4.75e+05]
146
The five (5) quantile classes, ranging from 0–1,255, 1,255–7,121, 7,121–18,488, 18,488–39,757, and 39,757–475,055 trips, confirm a highly right-skewed distribution where few hexagons capture exceptionally high demand while most remain low. The weekday morning peak exhibits a polarized pattern with many cells in both the lowest and highest classes, reflecting concentrated commuter flows toward major employment corridors. The afternoon peak displays a flatter distribution, indicating a more spatially dispersed pattern of return trips across the network. The weekend and holiday period demonstrates the most balanced spread, showing widespread moderate activity consistent with leisure-oriented mobility. Minor deviations from equal shares across classes arise from pooling but ensure consistent visual comparability across maps. Overall, morning demand is spatially focused, evening demand is diffused, and weekend demand is evenly distributed, establishing clear expectations for subsequent spatial statistical analysis of clustering and emerging hot spot behavior.
6.3 Peak-period visualisation
Each mandated peak-period map is produced as an independent, publication-quality figure using consistent cartographic settings. The national mainland boundary provides spatial reference, while analytical hexagons are shaded by total origin trips from the defined start-hour periods. Cell borders are omitted for visual clarity, and external legends enhance readability. A uniform sequential palette and the shared class intervals from Section 6.2 ensure identical colour meanings across all maps. Rendering each map separately prevents overlap and scaling distortions, maintaining analytical precision and reproducibility. The subsequent subsections present the code for the weekday morning, weekday afternoon, and weekend or holiday peaks, all based on the validated mainland hexagon geometry.
# --- Cartographic Style for Peak-Period Maps (AM, PM, Weekend/Holiday) -----# Shared layout template for consistency presentation# Function to generate a map with centered bold title and internal legendplot_peak_map <-function(data_sf, title_text, legend_title) {tmap_mode("plot") # static output for reproducible figuretm_shape(sg_outline) +tm_borders(col ="grey40", lwd =0.7) +# national boundary for spatial contexttm_shape(data_sf) +tm_polygons(col ="TRIPS",palette ="YlOrRd", # sequential colour schemebreaks = breaks_common, # shared 5-class quantilesborder.col ="grey", # light white edge for claritylwd =0.1,title = legend_title # legend title per map ) +tm_layout(main.title = title_text, # map title textmain.title.position ="center", # centre the titlemain.title.size =1.3, # slightly larger fontmain.title.fontface ="bold", # bold stylelegend.position =c("right", "bottom"), # legend inside map framelegend.bg.color ="white", # white background for readabilitylegend.frame =FALSE, # draw legend framelegend.text.size =0.65, # readable fontlegend.title.size =0.7,frame =TRUE, # draw neatline around mapouter.margins =0, # remove external paddinginner.margins =c(0.02, 0.02, 0.08, 0.02) # room for legend )}# --- Prepare sf objects for mapping ----------------------------------------# Attach geometry for each period total so that they can be plotted by tmap.# Weekday morning peak (start-hours 06 07 08)map_am <- hexagon_active %>%# take validated mainland hexagons dplyr::select(HEX_ID, geometry) %>%# keep only id and geometry dplyr::left_join(tp_am, by ="HEX_ID") %>%# attach morning totals sf::st_as_sf() # ensure sf class for tmap# Weekday afternoon peak (start-hours 17 18 19)map_pm <- hexagon_active %>% dplyr::select(HEX_ID, geometry) %>% dplyr::left_join(tp_pm, by ="HEX_ID") %>% sf::st_as_sf()# Weekend and holiday peak (start-hours 11 – 19)map_wk <- hexagon_active %>% dplyr::select(HEX_ID, geometry) %>% dplyr::left_join(tp_wk, by ="HEX_ID") %>% sf::st_as_sf()# --- Generate maps ----------------------------------------------------------# 1) Weekday Morning Peakplot_peak_map( map_am,title_text ="Weekday Morning Peak Origin Trips (Start-hours 06, 07, 08)",legend_title ="Trips (AM Peak)")
# 3) Weekend and Holiday Peakplot_peak_map( map_wk,title_text ="Weekend and Holiday Peak Origin Trips (Start-hours 11–19)",legend_title ="Trips (Weekend/Holiday Peak)")
The above three (3) maps show distinct spatial and temporal variations in passenger trip origins across Singapore. During weekday morning peaks (06–08), the highest concentrations of trips appear in residential heartlands such as Woodlands, Yishun, Sengkang, Tampines, Bedok, Bukit Panjang, and Jurong West, reflecting commuters travelling from home areas toward employment centres. In contrast, the weekday afternoon peaks (17–19) display a reversed pattern with intensified activity around the Central Region, notably along Orchard, Downtown Core, and Kallang, corresponding to return trips from workplaces to suburban residences. On weekends and holidays (11–19), the trip origins are more spatially dispersed but remain prominent in retail and leisure corridors such as Orchard, Marina Bay, and East Coast, as well as regional town centres. The consistent clustering along major transport corridors confirms strong home-to-work commuting patterns on weekdays and leisure-driven movement on weekends, illustrating the temporal shift in public transport demand across Singapore’s mainland.
7 Spatial Weights Construction
Spatial weights matrices define the strength of spatial relationships between neighbouring analytical units. In this study, spatial weights were constructed for the 375 m hexagonal tessellation to capture inter-hexagonal interactions in bus trip intensity. These weights quantify how one hexagon’s attribute value (for instance, total boardings) relates to those of surrounding hexagons. Following the methodology in Hands-on Exercise 4, 2 structures were generated — Queen contiguity and distance-band neighbours — each reflecting a different spatial interaction assumption. Contiguity weights consider only physical adjacency (shared borders or corners), whereas distance-band weights consider all cells within a specified metric distance. This hybrid design allows subsequent spatial autocorrelation analyses to be sensitive both to geographical connectivity and to spatial proximity. All computations employed spdep, sf, and tmap, with geometries validated and projected in SVY21 (Singapore Meters) to ensure Euclidean accuracy. The final selected structure, a 750 m distance-band model, provides complete spatial connectivity while preserving realistic local influence.
7.1 Queen Contiguity Weights
This section establishes the spatial contiguity framework required for subsequent spatial autocorrelation analysis in this study. The code chunk first retains only the unique hexagon identifier and geometry from the active analytical grid. Using poly2nb(), a Queen contiguity structure is generated, recognising neighbours that share either edges or corners. The code then counts and reports any zero-degree (isolated) hexagons to ensure there are no disconnected spatial units that could bias later calculations. Next, the neighbour list is converted into a row-standardised spatial-weights object (nb2listw) under the “W” style, allowing inclusion of cells with no neighbours (zero.policy = TRUE). Finally, the summary() function provides a diagnostic summary of neighbour counts for the Queen configuration, confirming the structural integrity of the spatial-weights matrix used in subsequent analyses.
# --- 7.2 Queen contiguity ---------------------------------------# retain only identifier and geometryhex <- hexagon_active |> dplyr::select(HEX_ID, geometry)# compute queen (edge or corner) neighboursnb_q <- spdep::poly2nb(as_Spatial(hex), queen =TRUE)# count zero-degree (isolated) cells for diagnosticsn_zero_q <-sum(card(nb_q) ==0) cat("\n Queen:", n_zero_q, "isolated")
# Inspect the numeric of Queen contiguitysummary(lengths(nb_q))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 4.000 5.000 4.793 6.000 6.000
# Inspect the characteristics of row-standardised weightlw_q_W
Characteristics of weights list object:
Neighbour list object:
Number of regions: 826
Number of nonzero links: 3958
Percentage nonzero weights: 0.5801171
Average number of links: 4.791768
1 region with no links:
826
6 disjoint connected subgraphs
Weights style: W
Weights constants summary:
n nn S0 S1 S2
W 825 680625 825 384.3206 3335.117
The results confirm that the Queen contiguity network was successfully generated with 826 hexagonal regions, of which one cell is isolated and lacks neighbours. The network comprises 3958 nonzero spatial links, representing about 58% of all possible neighbour connections. Each hexagon has an average of 4.79 neighbours, ranging from 1 to 6, consistent with the expected geometry of a well-structured hexagonal grid. The presence of six disjoint connected subgraphs suggests minor fragmentation, likely at boundary or offshore areas. After applying row standardisation (style = “W”), the resulting spatial-weights matrix retains the same neighbour relationships while normalising each row to ensure equal influence across all hexagons. The output summary of weights constants further validates the integrity and completeness of the Queen-based spatial-weight structure, which is ready for subsequent spatial autocorrelation analysis.
7.2 Distance-Band Weights
While contiguity captures geometric adjacency, it can fragment edge cells and isolate coastal zones. Therefore, we extend to a distance-band weights matrix, where hexagons within a specified metric threshold are considered neighbours. Candidate distance thresholds were derived empirically from the first-nearest-neighbour distance distribution of centroid-to-centroid separations. Multiple quantiles (0.75–0.99) were evaluated using spdep::dnearneigh(), and the smallest upper bound yielding a single connected component was selected. For this study, that distance was 750 meters — approximately twice the hexagon diameter—ensuring topological completeness without over-connecting remote cells. This data-driven method balances the need for analytic connectivity and local spatial realism.
7.2.1 Derive first-NN distances
Creates a working hexagon layer (HEX_ID, geometry) and takes one representative point per hex (st_point_on_surface). Computes the full pairwise distance matrix D, removes units, sets the diagonal to Inf, then extracts each hexagon’s nearest-neighbour distance as d1 = apply(D, 1, min). Non-positive or non-finite values are dropped and a summary of d1 is shown. This vector is the empirical basis for all subsequent distance-band choices.
Min. 1st Qu. Median Mean 3rd Qu. Max.
750.0 750.0 750.0 751.8 750.0 2250.0
The summary output shows that the nearest-neighbour (first-NN) distances among hexagons mostly converge at 750 m, with the minimum, first quartile, median, and third quartile all equal to this value. The mean is slightly higher at ≈751.8 m, indicating minimal variation, while the maximum distance of 2250 m reflects a few isolated or edge hexagons that are farther from their nearest neighbours. Overall, the results confirm that the hexagonal grid is uniformly spaced, producing consistent neighbour proximity across most of Singapore’s analytical surface. The isolated large distances likely occur at boundary areas where surrounding hexagons are fewer or disconnected. This outcome validates that 750 m is the typical inter-centroid spacing within the active hexagon system, making it a logical candidate for defining the initial distance threshold in subsequent neighbour and spatial-weight construction steps.
7.2.2 Scan candidate thresholds
Builds candidate upper-band distances from selected quantiles of d1 (75% to 99%). Converts points to numeric coordinates and defines pick_upper(ub) which forms dnearneigh(d1=0, d2=ub) and records the number of connected components. Applying it over all candidates yields scan_df, a table of (upper, components) for diagnosing when the graph becomes connected.
The table shows the number of connected components across different candidate upper-band thresholds, derived from the 75th to 99th percentiles of first-neighbour distances. All thresholds below the 99th percentile produce 79 disconnected components, meaning most hexagons remain isolated at shorter distances. When the upper limit reaches the 99th percentile (750 m), the number of components drops sharply to 6, indicating that larger connection bands begin merging separate clusters into broader connected regions. However, the network is still not fully unified into a single component. This outcome suggests that a threshold slightly above 750 m may be required to achieve complete connectivity across all hexagons. The result helps determine the minimal distance needed for every analytical cell to have at least one valid neighbour, ensuring that subsequent spatial-weight matrices are fully connected for reliable spatial dependence testing.
7.2.3 Select minimal connected threshold
Chooses the smallest upper-band that yields a single connected component (components == 1). If such a row exists, ubest is its upper. If none exists, it safely falls back to the largest tested upper-band. The result ubest is the data-driven distance cut-off to ensure a connected neighbour graph with minimal stretching.
The output identifies 750 meters as the optimal upper distance threshold (ubest) that ensures all hexagons form a single connected spatial network. This means that when hexagon centroids within 750 m of each other are treated as neighbours, every analytical unit becomes part of a continuous contiguity structure with no isolated subgraphs. The conditional check confirms that no smaller distance achieves full connectivity, so 750 m is the minimal connected threshold. Selecting this value provides a balance between spatial precision and network completeness, avoiding unnecessary link expansion while maintaining analytical integrity. This threshold will be used in the next step to construct the final distance-based spatial weight matrix, ensuring consistent neighbour relationships for spatial autocorrelation analysis.
7.2.4 Build neighbour list and weights
Constructs the final distance-band neighbours with dnearneigh(d1=0, d2=ubest). Summaries of neighbour counts are printed. Two weighting schemes are prepared: (i) binary weights (style=“B”, zero.policy=TRUE) via nb2listw; and (ii) inverse-distance weights by computing nbdist on the coordinates, inverting distances (\(1/v\)), and passing them as glist to nb2listw with style="B" and zero.policy=TRUE.
# --- 7.3.4 Build neighbour list and weights -------------------------------nb_d <- spdep::dnearneigh(coords_final, d1 =0, d2 = ubest)summary(lengths(nb_d))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 4.000 5.000 4.793 6.000 6.000
The summary output shows that each hexagon has between 1 and 6 neighbours, with an average of about 4.79, indicating a well-connected grid structure. The warning — “neighbour object has 6 sub-graphs” — signals that the network remains divided into six disconnected spatial clusters, meaning not all hexagons are part of a single continuous graph. This segmentation is expected in edge areas or offshore regions where some cells lack adjacent counterparts. Two spatial-weight schemes were successfully created: binary weights, where all neighbours are equally weighted, and inverse-distance weights, where closer neighbours exert stronger influence. Together, these provide flexible foundations for spatial autocorrelation and clustering analyses in later steps. Despite the warning, the results confirm the robustness of the 750 m threshold, ensuring most hexagons have valid neighbour relationships while maintaining high spatial coherence across Singapore’s analytical grid.
7.3 Connectivity diagnostics and visualisation
In this section, we will define a custom function nb_to_lines() to convert neighbour relationships from the spatial-weights object into line geometries connecting each pair of neighbouring hexagons. It loops through each neighbour pair, constructs LINESTRING geometries between their centroid coordinates, and compiles them into a simple feature (sf) object named edges_db. The map visualisation then overlays three layers: Singapore’s outline (sg_outline), the analytical hexagons (hex_proj), and the neighbour links (edges_db) drawn in dark orange to represent connectivity. Additional map elements include a title displaying the distance-band threshold (ubest), a legend explaining “Neighbour link,” and customised styles for borders, line width, and layout positioning. Together, the code provides a clear visual representation of spatial relationships among hexagons, confirming how the distance-based network connects neighbouring cells across the study area.
The iterative map visualises the distance-band neighbour network for Singapore’s analytical hexagons using an upper threshold of 750 meters. Each orange line represents a valid neighbour connection between adjacent hexagons within this distance, forming a continuous lattice across most of the mainland. The dense network in central and eastern regions confirms strong spatial connectivity, while small gaps or missing links near coastlines and peripheral zones indicate edge effects where fewer neighbouring cells exist. Overall, the 750 m threshold effectively balances network coverage and spatial precision—ensuring nearly all hexagons are linked without excessive overlap. This visual validation confirms that the spatial-weight structure is coherent and ready for subsequent spatial autocorrelation and clustering analyses.
8 Global Measure of Spatial Autocorrelation Analysis
This section investigates the global spatial autocorrelation of origin trip intensity across the analytical hexagon grid. Rather than simply examining how values vary geographically, the analysis quantifies the extent to which neighbouring cells exhibit similar or dissimilar trip counts — that is, whether the observed pattern represents clustering, randomness, or dispersion. The global Moran’s I statistic is used as the principal measure to capture this spatial dependency, providing a single index summarising overall correlation across the study region. A positive Moran’s I indicates that cells with high (or low) trip values are located near others of similar magnitude, implying spatial clustering. Conversely, a negative Moran’s I denotes spatial dispersion where high and low values alternate spatially. The test operates on the row-standardised Queen contiguity weights constructed earlier, ensuring that each cell’s contribution is proportional to its number of neighbours. The resulting Moran’s I value, standard deviate, and \(p\)-value are derived under a randomisation hypothesis, allowing formal inference against the null of spatial independence. Through this procedure, we obtain a statistically robust understanding of whether trip distributions during peak periods are spatially structured, thereby establishing the foundation for subsequent LISA and hot/cold-spot analyses.
8.1 Global Moran’s I statistical test for all peak periods
This section performs the Global Moran’s I test independently for the three defined peak periods: weekday morning, weekday afternoon, and weekend/holiday. Using the previously derived Queen-contiguity spatial-weights matrix, each hexagonal layer is tested under the null hypothesis of spatial randomness. The moran.test() computes the observed Moran’s I, expected value, variance, \(z\)-score, and corresponding \(p\)-value. A positive and significant Moran’s I (\(p < 0.001\)) confirms that the trip intensities are spatially clustered during that period, while an insignificant result suggests no global pattern. Comparing Moran’s I values across all periods reveals whether clustering strength is consistent or time-specific, providing evidence of temporal variation in spatial dependency.
# --- 8.1 Global Moran’s I (AM, PM, Weekend) ---------------------------------# helper: make the numeric vector aligned to the weights (reuses Section 6/7 objects)make_x <-function(trip_table) { hexagon_active %>% dplyr::select(HEX_ID, geometry) %>% dplyr::left_join(trip_table %>% dplyr::select(HEX_ID, TRIPS), by ="HEX_ID") %>% dplyr::mutate(Trips = dplyr::coalesce(TRIPS, 0)) %>% {.$Trips}}x_am <-make_x(tp_am)x_pm <-make_x(tp_pm)x_wk <-make_x(tp_wk)# Moran's I test under randomisationmt_am <- spdep::moran.test(x_am, listw = lw_d_B, zero.policy =TRUE, na.action = na.omit)mt_pm <- spdep::moran.test(x_pm, listw = lw_d_B, zero.policy =TRUE, na.action = na.omit)mt_wk <- spdep::moran.test(x_wk, listw = lw_d_B, zero.policy =TRUE, na.action = na.omit)mt_am; mt_pm; mt_wk # print results
Moran I test under randomisation
data: x_am
weights: lw_d_B
n reduced by no-neighbour observations
Moran I statistic standard deviate = 16.367, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.3627974234 -0.0012135922 0.0004946119
Moran I test under randomisation
data: x_pm
weights: lw_d_B
n reduced by no-neighbour observations
Moran I statistic standard deviate = 7.1164, p-value = 5.538e-13
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.1542609092 -0.0012135922 0.0004773017
Moran I test under randomisation
data: x_wk
weights: lw_d_B
n reduced by no-neighbour observations
Moran I statistic standard deviate = 10.997, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.2416785594 -0.0012135922 0.0004878754
The randomisation tests show strong positive spatial autocorrelation for trip counts in all periods using the lw_d_B weights. For weekday morning (06–09), Moran’s I = 0.3628, \(z\) = 16.367, \(p\) < 2.2×10⁻¹⁶, decisively rejecting spatial randomness. For weekday afternoon (17–20), Moran’s I = 0.1543, \(z\) = 7.116, \(p\) = 5.54×10⁻¹³, indicating weaker—but still highly significant—clustering. For weekend/holiday (11–19), Moran’s I = 0.2417, \(z\) = 10.997, \(p\) < 2.2×10⁻¹⁶, implying moderate clustering. The “Expectation” of about −0.0012 is the theoretical mean of Moran’s I under the null (approximately zero for large samples), and the listed variances quantify sampling variability under randomisation. Together, large positive \(I\) values, very high \(z\)-scores, and extreme \(p\)-values confirm that neighbouring hexagons tend to have similar trip intensities, with clustering strongest in the morning, moderate on weekends, and weakest in the weekday afternoon.
8.2 Computing Monte Carlo Moran’s I
The code chunk below performs permutation test for Moran’s I statistic by using moran.mc() of spdep. A total of 999 simulation will be performed.
Monte-Carlo simulation of Moran I
data: x_am
weights: lw_d_B
number of simulations + 1: 1000
statistic = 0.3628, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
Monte-Carlo simulation of Moran I
data: x_pm
weights: lw_d_B
number of simulations + 1: 1000
statistic = 0.15426, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
Monte-Carlo simulation of Moran I
data: x_wk
weights: lw_d_B
number of simulations + 1: 1000
statistic = 0.24168, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
The Monte-Carlo simulation results confirm significant positive spatial autocorrelation in trip intensities across all time periods, using 999 random permutations. For the weekday morning (06–09), Moran’s I = 0.3628 with p = 0.001, indicating strong clustering of high-trip hexagons—typical of concentrated commuting activity. During the weekday afternoon (17–20), Moran’s I decreases to 0.1543, suggesting weaker but still significant clustering as trip origins become more spatially dispersed. On weekends and holidays (11–20), Moran’s I = 0.2417, showing moderate clustering, consistent with concentrated leisure or retail-related travel zones. The consistently low \(p\)-values (< 0.001) across all three tests reject spatial randomness, confirming that trip distributions remain spatially dependent, though clustering strength varies with the time of day.
8.3 Visualising Monte Carlo Moran’s I
It is always a good practice for us the examine the simulated Moran’s I test statistics in greater detail. This can be achieved by plotting the distribution of the statistical values as a histogram by using the code chunk below.
The permutation histograms illustrate the simulated distributions of Moran’s I under spatial randomness, with the red line marking the observed Moran’s I for each time period. In all three plots, the observed values lie far to the right of their respective simulated distributions, indicating strong positive spatial autocorrelation. The weekday morning (06–09) shows the largest deviation, confirming intense clustering of trip origins during peak commute hours. The weekday afternoon (17–20) displays a smaller separation, reflecting weaker clustering and greater dispersion of trips. The weekend/holiday (11–20) pattern falls between the two, showing moderate clustering likely related to leisure activity centres. Since the observed Moran’s I values are far outside the random distribution in every case, spatial randomness is rejected with high confidence, validating significant spatial dependence in all time periods.
8.4 Global Moran’s I spatial correlogram
This section constructs the spatial correlogram of Moran’s I to assess how spatial autocorrelation changes with increasing distance. The correlogram plots Moran’s I values against successive distance bands derived from the spatial-weights structure, effectively tracing how similarity between cells weakens or reverses as neighbourhoods expand. A steep decline of Moran’s I toward zero implies a rapid loss of spatial dependency, whereas sustained positive values across several lags suggest broader clustering effects. This diagnostic complements the global Moran’s I test by revealing the scale and range of spatial association, thus offering deeper insight into the spatial processes shaping trip-generation patterns across different peak periods.
# --- 8.2 Spatial correlogram (Moran’s I by lag order) -----------------------library(spdep)# Use the neighbour list from Section 7 (nb_d) and W-style like in classMIcorr_am <- spdep::sp.correlogram(nb_d, x_am, order =6, method ="I", style ="W", zero.policy =TRUE)MIcorr_pm <- spdep::sp.correlogram(nb_d, x_pm, order =6, method ="I", style ="W", zero.policy =TRUE)MIcorr_wk <- spdep::sp.correlogram(nb_d, x_wk, order =6, method ="I", style ="W", zero.policy =TRUE)# plot each correlogram and print the tables (as shown in the notes)plot(MIcorr_am, main ="Moran’s I correlogram – AM 06–09"); print(MIcorr_am)
The spatial correlograms and tables show how Moran’s I changes across 6 lags for the 3 periods. In every case, lag 1 is strongly positive and highly significant, meaning neighbouring hexagons (at the closest lag) have similar trip intensities. The weekday morning (06–09) exhibits the largest short-range dependence (\(I\) ≈ 0.343, \(z\) ≈ 14.6, \(p\) < 2.2e-16). Its autocorrelation declines rapidly with distance, becoming weak by lags 4–6; significance disappears from lag 5 onward, indicating a short spatial reach around the morning hotspots. The weekday afternoon (17–20) starts lower (\(I\) ≈ 0.152, \(z\) ≈ 6.64, \(p\) ≈ 3.1e-11) and decays steadily; lags 2–4 remain positive and significant but smaller, lag 5 is not significant, and lag 6 shows a small yet significant positive value, implying faint longer-range structure. Weekend/holiday (11–20) sits between morning and afternoon: strong at lag 1 (\(I\) ≈ 0.230, \(z\) ≈ 9.88, \(p\) < 2.2e-16), still significant through lag 4, not significant at lag 5, and again modestly positive at lag 6. Overall, spatial dependence is strongest and most localized in the morning, weakest in the afternoon, and moderate on weekends. The attenuation patterns confirm that clustering diminishes with distance, with occasional long-range remnants at the farthest lag.
9 Local Measure of Spatial Autocorrelation Analysis
9.1 Computing Local Moran’s I
In this section, a custom function is designed to calculate the Local Moran’s I statistic for each analytical hexagon. The function integrates trip intensity data with spatial geometry, replaces missing values with zeros, and performs 999 random permutations using distance-based spatial weights. This computation produces multiple diagnostic outputs, including local Moran’s I values, variances, \(z\)-scores, and \(p\)-values. Separate analyses are then performed for morning, afternoon, and weekend periods to capture time-based spatial clustering variations.
nsim <-999# keep same as earlier sections# helper: join trips to hexagons, then compute local Moran's I (sfdep)compute_lisa_sf <-function(trip_table) { var_df <- trip_table %>%select(HEX_ID, TRIPS) %>%rename(Trips = TRIPS) lisa <- hexagon_active %>%select(HEX_ID, geometry) %>%left_join(var_df, by ="HEX_ID") %>%mutate(Trips =as.numeric(coalesce(Trips, 0))) %>%# numeric + fill edges with 0mutate(local_moran =local_moran(x = Trips,nb = nb_d, # nb list (Section 7)wt = lw_d_B$weights, # <-- pass the WEIGHTS LIST, not the listwnsim = nsim,zero.policy =TRUE ),.before =1 ) %>%unnest(local_moran) # adds ii, e_ii, var_ii, z_ii, p_ii lisa}# build per period (reusing tp_* from Section 6)lisa_am <-compute_lisa_sf(tp_am)lisa_pm <-compute_lisa_sf(tp_pm)lisa_wk <-compute_lisa_sf(tp_wk)# Inspect the data structure glimpse(lisa_am)
The printed data frames are the per-hexagon outputs from the Local Moran’s I calculation (826 rows). For each hexagon, ii is the observed local Moran’s I, e_ii is its null expectation (near −0.0012), and var_ii the randomisation variance. The z_ii column standardises ii; large positive values indicate clustering of similar high or low trip counts, while negative values indicate spatial dispersion or outliers. Columns p_ii (analytic), p_ii_sim (permutation), and p_folded_sim report significance; many entries are extremely small (often ≤0.01), signalling numerous statistically significant local associations. The categorical labels mean, median, and pysal assign each hexagon to a LISA quadrant; in this excerpt, labels are overwhelmingly Low–Low, meaning hexagons with low trips surrounded by neighbours that are also low. Diagnostic fields skewness and kurtosis describe the shape of the permutation distribution for each unit. Identifiers HEX_ID, Trips, and geometry link results back to the grid and counts for mapping. Overall, the tables show widespread significant Low–Low clustering with pockets of other patterns flagged by positive z-scores and tiny \(p\)-values.
9.1.1 Mapping the Local Moran’s I
This part focuses on visualising the computed Local Moran’s I results as thematic maps. A mapping function is defined to display the spatial variation of Moran’s I values across Singapore, using a continuous colour scale to represent clustering strength. The maps are uniformly formatted with consistent legends and titles, allowing for direct comparison between time periods and facilitating the detection of spatial heterogeneity within the analytical grid.
The Local Moran’s I maps provide spatially explicit insights into the degree of clustering or dispersion of trip origins during different time periods. In the weekday morning (06–09), pronounced positive clustering (dark blue hexagons) appears in the western, northeastern, and parts of the central regions, suggesting areas where high trip intensities are surrounded by similarly high neighbours. These patterns correspond to residential zones with strong outbound commuting flows. Negative values (orange tones) indicate isolated cells or transitional areas where trip activity diverges from nearby cells, suggesting spatial discontinuities likely caused by mixed land use or low transport demand. During the weekday afternoon (17–20), clustering weakens, with fewer high-positive Moran’s I values, indicating that evening trips are more spatially dispersed as commuters return to various home locations. In the weekend/holiday (11–19) period, moderate positive clusters persist mainly in the southern and central regions, reflecting concentrated leisure or shopping-related trips in activity hubs like Orchard or Marina areas. Overall, spatial autocorrelation is strongest during weekday mornings and declines across later periods, showing that Singapore’s trip-generation patterns transition from structured, residential-based flows to more scattered, destination-driven movements. This implies the morning peak remains the most spatially organised travel window, essential for targeting transport interventions and infrastructure optimisation.
9.1.2 Mapping Local Moran’s I \(p\)-values
Here, the emphasis shifts to mapping the statistical significance of local spatial relationships. \(P\)-values derived from Local Moran’s I computations are visualised using colour gradients, where darker tones highlight stronger evidence of spatial dependence. These maps reveal which areas exhibit meaningful local clustering and which remain random. Comparing across time periods provides insights into how the significance of clustering changes during different travel intervals throughout the day.
tmap_mode("plot")map_localP <-function(sf_obj, title_text) {tm_shape(sg_outline) +tm_borders(col ="grey40", lwd =0.7) +tm_shape(sf_obj) +tm_fill(col ="p_ii",fill.scale =tm_scale_intervals(breaks =c(-Inf, 0.001, 0.01, 0.05, 0.10, Inf),values ="brewer.Reds" ),fill.legend =tm_legend(title ="p-value"),showNA =FALSE, # hide "Missing" in legendcolorNA =NA# transparent for NA hexes ) +tm_borders(fill_alpha =0.5, col ="grey80", lwd =0.4) +tm_layout(main.title = title_text,main.title.position ="center",main.title.fontface ="bold",legend.position =c("right", "bottom"),legend.frame =TRUE )}# Run all 3 mapsmp_am <-map_localP(lisa_am, "p-values of Local Moran's I — Weekday Morning (06–09)")mp_pm <-map_localP(lisa_pm, "p-values of Local Moran's I — Weekday Afternoon (17–20)")mp_wk <-map_localP(lisa_wk, "p-values of Local Moran's I — Weekend/Holiday (11–20)")mp_am; mp_pm; mp_wk
The \(p\)-value maps visualise the statistical significance of local clustering across Singapore for different periods. In the weekday morning (06–09), numerous light-blue hexagons (p ≤ 0.2) dominate the map, especially in western, northeastern, and southeastern regions, signifying that local Moran’s I values are statistically significant and that spatial clustering of trip origins is unlikely to occur by chance. These significant cells correspond to areas with concentrated commuting activity, validating the earlier high Moran’s I findings. During the weekday afternoon (17–20), significant regions decrease slightly, reflected by more mid-blue shades (\(p\) ≈ 0.4–0.6), suggesting reduced clustering as trip destinations diversify across urban zones. The weekend/holiday (11–20) period shows the fewest highly significant cells, particularly in residential and industrial zones, where darker blue (\(p\) > 0.6) indicates random or weak spatial association. This pattern aligns with more dispersed, leisure-based travel during non-work periods. Overall, the intensity and spatial coverage of low \(p\)-values decline from morning to weekend, implying that weekday peak-hour movements are spatially structured, while off-peak and weekend activities exhibit more randomness. The outcome underscores that spatial dependence in trip generation is temporally dynamic—strongest during structured commuting peaks and weakest during flexible leisure travel—informing transport planners about when and where spatial targeting of travel interventions is most relevant.
9.1.3 Mapping both Local Moran’s I values and \(p\)-values
This segment combines two key perspectives by displaying Moran’s I values and their associated \(p\)-values side by side. The paired visualisation allows simultaneous interpretation of spatial clustering strength and statistical reliability. By comparing both maps together, analysts can easily distinguish significant clusters from random fluctuations. The approach improves interpretability and ensures that subsequent discussions are grounded in statistically validated spatial patterns across all travel periods.
tmap_arrange(map_localI(lisa_am, "Local I — Morning (06–09)"),map_localP(lisa_am, "p-value — Morning (06–09)"),ncol =2, asp =1)
tmap_arrange(map_localI(lisa_pm, "Local I — Afternoon (17–20)"),map_localP(lisa_pm, "p-value — Afternoon (17–20)"),ncol =2, asp =1)
tmap_arrange(map_localI(lisa_wk, "Local I — Weekend/Holiday (11–20)"),map_localP(lisa_wk, "p-value — Weekend/Holiday (11–20)"),ncol =2, asp =1)
The combined Local Moran’s I and \(p\)-value maps reveal both the strength and statistical reliability of spatial clustering across Singapore’s trip-origin data. During the weekday morning (06–09), areas with high positive Local I (dark blue on the left) coincide with light-blue \(p\)-values (right map), especially in the west and northeast. These zones represent statistically significant clusters (p < 0.05) where strong spatial dependence exists—typically high-trip cells surrounded by high-trip neighbours. In contrast, orange cells (negative I) appear near the central fringe, indicating dissimilar neighbourhood effects such as isolated business or institutional zones. In the weekday afternoon (17–20), overall clustering intensity weakens, and significant cells shrink mainly to residential–commercial transition belts, reflecting dispersed evening travel as commuters return home. By weekend/holiday (11–20), high I values persist only around southern activity hubs, but most regions exhibit non-significant p-values (dark blue on the right), confirming randomised spatial patterns during leisure hours. The progressive reduction of significant low-p clusters from morning to weekend implies that trip generation is most spatially structured during work-related commuting peaks but becomes increasingly stochastic when travel is discretionary. For transport planners, this pattern highlights when and where spatial dependencies dominate mobility flows—vital for prioritising public-transport capacity, infrastructure optimisation, and land-use coordination during critical morning peaks.
9.2 Preparing and visualising LISA Map
At this stage, the analysis transitions toward visual diagnostics by preparing data for Moran scatterplots. The process involves computing spatially lagged trip counts, ensuring new variables are properly reassigned, and plotting the relationships between each hexagon’s trip value and its neighbours. The resulting scatterplots highlight the overall correlation structure and clustering direction, helping confirm that local spatial dependencies identified earlier are both consistent and statistically meaningful across all study periods.
9.2.1 Plotting Moran Scatterplot
Moving from maps to diagnostics, this section constructs Moran scatterplots for each time period. A helper adds the spatial lag of Trips directly from the existing listw, ensuring consistency with the distance band weights. The plots place Trips on the x axis against their spatial lags on the y axis and fit an OLS line to reveal spatial association. Points display individual hexagons, while the fitted slope visualises positive or negative dependence prior to standardisation.
### Add spatial lag of Trips for scatterplots (no unnesting, no wrapper)add_lag <-function(sf_obj) {# compute lag directly from the listw object you already built sf_obj$lag_Trips <-as.numeric( spdep::lag.listw(lw_d_B, sf_obj$Trips, zero.policy =TRUE, NAOK =TRUE) ) sf_obj}# REASSIGN so the new column stickslisa_am <-add_lag(lisa_am)lisa_pm <-add_lag(lisa_pm)lisa_wk <-add_lag(lisa_wk)moran_scatter <-function(sf_obj, title_text) {ggplot(sf_obj, aes(x = Trips, y = lag_Trips)) +geom_point() +geom_smooth(method ="lm", se =FALSE, color ="red") +labs(x ="Trips", y ="Spatial Lag of Trips", title = title_text) +theme_minimal()}ms_am <-moran_scatter(lisa_am, "Moran Scatterplot — Weekday Morning (06–09)")ms_pm <-moran_scatter(lisa_pm, "Moran Scatterplot — Weekday Afternoon (17–20)")ms_wk <-moran_scatter(lisa_wk, "Moran Scatterplot — Weekend/Holiday (11–20)")ms_am; ms_pm; ms_wk
he Moran scatterplots illustrate the spatial relationship between trip intensity (x-axis) and the spatial lag of trips (y-axis), showing how local trip counts correlate with neighbouring areas across three time periods. The weekday morning (06–09) plot exhibits a clear upward trend with a steep red regression line, indicating strong positive spatial autocorrelation. This means hexagons with high trip numbers are surrounded by similarly high-trip neighbours, reflecting concentrated commuting activity in key employment and transport nodes. During the weekday afternoon (17–20), the slope remains positive but flatter, suggesting a weaker spatial dependency. This reflects more spatial dispersion as commuters return home, reducing the dominance of clustered high-trip zones. On weekends/holidays (11–20), the scatter pattern stays positively sloped but becomes even more diffuse, indicating lower clustering and greater variability in trip patterns, typical of leisure or non-routine movements. The persistence of a positive slope across all periods confirms that spatially proximate hexagons tend to share similar trip magnitudes, but the declining gradient from morning to weekend highlights a shift from structured, work-related clustering to spatially scattered travel behaviour. These patterns underscore that spatial dependency in mobility is strongest during structured weekday peaks and weakest when travel choices are discretionary, providing critical insights for adaptive transport planning and congestion management.
9.2.2 Plotting Moran scatterplot with standardised variable
This section enhances interpretability by standardising both trip counts and their spatial lags before plotting. Each hexagon is classified into one of four categories—High–High, Low–Low, High–Low, or Low–High—according to \(z\)-scores. Distinct colours are applied to visualise these quadrants, making cluster and outlier patterns easily recognisable. The standardised Moran scatterplots provide an intuitive visual explanation of how different areas relate spatially, illustrating both concentration zones and transition boundaries across Singapore.
#### 10.4.4.2 Plotting Moran scatterplot with standardised variable ----# 1) Standardise Trips and its spatial lag (center/scale, return plain numeric)standardise_for_scatter <-function(sf_obj){ sf_obj |> dplyr::mutate(z_Trips =as.vector(scale(Trips)),z_lag_Trips =as.vector(scale(lag_Trips)) )}lisa_am <-standardise_for_scatter(lisa_am)lisa_pm <-standardise_for_scatter(lisa_pm)lisa_wk <-standardise_for_scatter(lisa_wk)# 2) (Guard) ensure LISA quadrant labels `mean` exist as in 10.4.4.1# If your earlier chunk already created `mean`, this block does nothing.ensure_quadrants <-function(sf_obj){if (!"mean"%in%names(sf_obj)) { x0 <-mean(sf_obj$Trips, na.rm =TRUE) y0 <-mean(sf_obj$lag_Trips, na.rm =TRUE) sf_obj <- dplyr::mutate( sf_obj,mean = dplyr::case_when( Trips >= x0 & lag_Trips >= y0 ~"High-High", Trips < x0 & lag_Trips < y0 ~"Low-Low", Trips < x0 & lag_Trips >= y0 ~"Low-High",TRUE~"High-Low" ),mean =factor(mean, levels =c("High-High","Low-Low","Low-High","High-Low")) ) } sf_obj}lisa_am <-ensure_quadrants(lisa_am)lisa_pm <-ensure_quadrants(lisa_pm)lisa_wk <-ensure_quadrants(lisa_wk)# 3) Standardised Moran scatter (points + black OLS; dashed mean lines)library(ggplot2)moran_scatter_std <-function(sf_obj, title_text){ggplot(sf_obj, aes(x = z_Trips, y = z_lag_Trips, color = mean)) +geom_point(size =2) +geom_smooth(method ="lm", se =FALSE, color ="black") +geom_hline(yintercept =mean(sf_obj$z_lag_Trips, na.rm =TRUE), lty =2) +geom_vline(xintercept =mean(sf_obj$z_Trips, na.rm =TRUE), lty =2) +scale_color_manual(values =c("High-High"="red","Low-Low"="blue","Low-High"="lightblue","High-Low"="pink" ) ) +labs(x ="Standardised Trips",y ="Standardised Spatial Lag of Trips",title = title_text ) +theme_minimal()}ms_am_z <-moran_scatter_std(lisa_am, "Standardised Moran Scatter — Weekday Morning (06–09)")ms_pm_z <-moran_scatter_std(lisa_pm, "Standardised Moran Scatter — Weekday Afternoon (17–20)")ms_wk_z <-moran_scatter_std(lisa_wk, "Standardised Moran Scatter — Weekend/Holiday (11–20)")ms_am_z; ms_pm_z; ms_wk_z
The standardised Moran scatterplots reveal the spatial association patterns of trip intensity by quadrant classification, distinguishing four spatial interaction types: High-High (HH), Low-Low (LL), High-Low (HL), and Low-High (LH). During the weekday morning (06–09), the majority of points fall in the High-High quadrant (red), confirming strong clustering of high-trip areas surrounded by similar high-trip neighbours — characteristic of concentrated commuting toward business and transport hubs. Low-Low zones (blue) dominate the opposite quadrant, indicating peripheral or residential regions with uniformly low trip activity. The weekday afternoon (17–20) plot still shows a positive spatial relationship, though the slope slightly flattens, suggesting more dispersed trip flows as outbound commuting begins. A few High-Low outliers (pink) appear, marking transition zones between dense business districts and adjacent low-activity cells. During the weekend/holiday (11–19), while High-High clusters remain visible, the overall distribution becomes broader with more scattered Low-High and High-Low points, reflecting irregular trip patterns driven by leisure or non-routine movement. The positive slope across all periods confirms persistent spatial autocorrelation, yet the reduction in clustering strength from weekday morning to weekend reflects a shift from structured work-related travel to more dispersed, flexible trip behaviour. This finding underscores the significance of weekday peaks for congestion management and highlights the spatial diffusion of mobility during non-working days, offering valuable input for dynamic transport planning and demand-based service allocation.
9.2.3 Preparing LISA Map Classes (\(p\) < 0.05)
Finally, the analysis establishes clear LISA map classifications based on significance thresholds. Each hexagon is assigned to a specific cluster category according to its Local Moran’s I value and \(p\)-value. Colours are systematically applied to differentiate strong clusters, weak outliers, and non-significant regions. The resulting maps visually summarise spatial autocorrelation outcomes for morning, afternoon, and weekend periods, helping to pinpoint areas of persistent high or low clustering within the study region.
The LISA Cluster Maps (\(p\) < 0.05) reveal clear temporal differences in local spatial autocorrelation across Singapore’s urban structure.
During **Weekday Morning (06–09), High–High (HH) clusters (dark red) are widespread and spatially concentrated across the central-northern, north-eastern, and western corridors. Prominent groupings appear around the northern central belt (likely Ang Mo Kio – Yishun region), another in the north-east zone (around Sengkang – Punggol), and a strong cluster in the western sector (Jurong – Clementi area). Smaller HH patches also appear in the south-central corridor, reflecting intensified morning mobility flows. The presence of Low–High (LH) cells (light blue) adjoining several HH patches suggests transitional areas—likely residential or mixed-use cells adjoining major trip origins or destinations. These spatial configurations indicate strong positive local autocorrelation, meaning high-trip hexagons are surrounded by similarly high neighbours, forming structured commuting corridors.
In Weekday Afternoon (17–20), HH clusters become sparser and more fragmented, persisting mainly in the north-east, west, and central-south belts. The contraction of HH zones and scattered LH outliers reflects spatial dispersion of evening travel, as movement decentralises from morning-dominated centres to more distributed return-trip destinations.
By Weekend/Holiday (11–20), HH clusters re-emerge more broadly across the northern, eastern, and south-central zones. These weekend clusters expand around previously active weekday corridors, implying recreational or shopping-related flows along regional centres. LH patches remain peripheral, delineating boundaries between active hubs and quieter suburbs.
10 Spatial Hot Spot and Cold Spot Analysis (Getis-Ord Gi*)
This section applies the Local Getis Ord Gi* statistic to identify statistically significant concentrations of high and low trip activities within the study area. The analysis builds upon the previously confirmed spatial weights and the hexagonal grid framework to evaluate whether the observed trip intensities form meaningful spatial clusters beyond what would be expected by random chance. Each analytical time window, including the weekday morning, weekday afternoon, and weekend or holiday periods, follows an identical spatial configuration to preserve methodological consistency. The workflow proceeds sequentially, beginning with the validation of spatial weight matrices, followed by the computation of local Gi* statistics and associated p values, and concluding with the visualization and classification of significant clusters. Every analytical step is executed reproducibly within the R environment through the sfdep and tmap packages, ensuring methodological alignment between statistical computation and geographic representation. The complete process transforms raw trip counts into interpretable spatial indicators of movement intensity, providing a rigorous foundation for understanding urban mobility patterns across temporal contexts.
10.1 Confirm weight objects (Include Self for Gi*)
The analysis begins by validating the spatial weight configurations required for the Gi* computation. The code first checks that the neighbour list object (nb_d_self) and its corresponding weight list (lw_d_B_self) exist and inherit from the correct nb and listw classes. The procedure explicitly includes each hexagon’s own location as a neighbour (include.self = TRUE) because Gi* requires self-inclusion to measure local concentration accurately. A sanity check confirms that the length of the neighbour structure matches the number of active hexagons (hexagon_active). These tests guarantee internal consistency between the spatial index and geometric dataset, preventing mismatched geometries or weight lengths that could invalidate later \(z\)-score and \(p\)-value calculations. By ensuring that the self-inclusive weights are both properly formatted and structurally aligned, this step forms the computational backbone for all subsequent spatial permutation tests and clustering assessments.
# Confirm 750 m weights from Section 7stopifnot(inherits(nb_d, "nb"),inherits(lw_d_B, "listw"))# Include-self versions (required for Gi*)nb_d_self <-if (exists("nb_d_self")) nb_d_self else spdep::include.self(nb_d)lw_d_B_self <-if (exists("lw_d_B_self")) lw_d_B_self else spdep::nb2listw(nb_d_self, style ="W", zero.policy =TRUE)
This section develops the computational tools required to implement the Local Getis Ord Gi* analysis reliably across multiple time periods. It introduces two essential helper functions that translate cleaned trip records into forms suitable for spatial computation. The first function ensures that all trip frequencies correspond exactly to their matching hexagonal cells within the analytical grid, producing a fully numeric vector that preserves spatial order and completeness. The second function performs the statistical calculation of the Local Gi* statistic using the previously validated spatial weights. Together, these procedures guarantee that each analytical dataset follows identical data structures and spatial relationships before entering the permutation process. The section also verifies that all outputs preserve consistent record counts and geometry alignment, ensuring no mismatch between the trip data and the spatial framework. Through these supporting functions, the analysis gains precision, reproducibility, and computational integrity, establishing a seamless bridge between data preparation and statistical computation in subsequent sections.
10.2.1 Convert trips to numeric vector aligned to grid
This section defines a helper routine that converts the tabulated trip data into a numeric vector perfectly aligned with the analytical hexagon grid. It performs a join between trip counts and geometry using the common hexagon identifier, replaces missing trip values with zero, and confirms that the resulting vector length matches the number of grid cells. Logical checks verify that all values are numeric and finite, ensuring that no data mismatch affects the spatial weights. The completed vector becomes the input for statistical analysis in later sections, providing a validated numerical representation of trip intensity that is consistent across all analytical windows.
trips_vec <-function(trip_table) { var_df <- trip_table %>%transmute(hex_id =tolower(HEX_ID),Trips =as.numeric(TRIPS)) x <- hexagon_active %>%transmute(hex_id =tolower(HEX_ID), geometry) %>%left_join(var_df, by ="hex_id") %>%mutate(Trips =replace_na(Trips, 0)) v <- x %>%st_drop_geometry() %>%pull(Trips)stopifnot(is.numeric(v),!any(is.na(v)),length(v) ==nrow(hexagon_active)) v}
10.2.2 Compute Local Gi* statistic
This section implements the computation of the Local Gi* statistic. The function compute_gistar uses the aligned trip vector and executes 999 random permutations with zero policy enabled to handle empty neighbour sets. The resulting \(z\)-score and two-sided \(p\) value quantify the degree and significance of clustering within each hexagon. The code then attaches the results to the analytical geometry, ensuring one record per spatial cell. Validation checks confirm identical row counts between statistical output and geometry, ensuring spatial consistency. This process yields an integrated spatial dataset containing both numeric and geometric information, ready for temporal comparison and visualisation.
After confirming the helper functions, the computation is applied to each temporal dataset representing the weekday morning, weekday afternoon, and weekend or holiday periods. The resulting outputs, labelled HCSA_am, HCSA_pm, and HCSA_wk, contain each hexagon’s Gi* \(z\)-score and simulated p value together with its temporal label. A diagnostic message confirms that all datasets contain identical numbers of observations and identical spatial structures. This repetition under uniform parameters allows precise comparison of spatial clustering strength between different time frames. The procedure captures temporal variation in the intensity of hot and cold clusters across Singapore, making it possible to observe whether strong morning concentrations persist into later periods or dissipate. The completion of this section produces three statistically comparable datasets that serve as the analytical input for the mapping process described in the following section.
This section transforms numerical results into spatially interpretable maps that visualise the Local Gi* \(z\)-scores. The mapping function uses tmap to display the Singapore outline, overlay the analytical hexagon grid, and fill each cell according to its \(z\)-score value. The chosen colour palette is diverging, with red tones representing positive clustering and blue tones representing negative clustering, while white represents values close to the mean. The map layout follows the standard formatting with bold titles, centred legends, and consistent margins to preserve visual uniformity across all time periods. Through these maps, users can immediately distinguish areas of elevated trip intensity from those with lower activity. The maps also provide a diagnostic check for spatial continuity, revealing whether the clustering patterns are geographically contiguous or dispersed. This step converts statistical output into an interpretable visual narrative that guides subsequent analysis of significance.
The three Local Gi* maps together present a detailed view of spatial clustering patterns of trip activity across different temporal contexts. During the weekday morning period, the map reveals intense hot spots represented by darker red hexagons located in the western, northeastern, and central regions. These areas correspond to concentrated morning departures associated with residential neighbourhoods where large numbers of commuters begin their daily journeys. The presence of light blue cells along the northern and southern corridors indicates localised cold spots, suggesting lower trip activity during this time window.
In the weekday afternoon period, several of the morning hot spots persist but appear more dispersed. This shift indicates that while core mobility zones remain active, their intensity diminishes as trip generation becomes more balanced between origins and destinations. The recurring pattern of both positive and negative \(z\)-scores suggests the redistribution of flows as commuters return from work.
On weekends and holidays, the pattern evolves further into more extensive red clusters in central and eastern regions, highlighting leisure-oriented and recreational travel rather than work-related mobility. The appearance of multiple moderate-intensity clusters across the island suggests more spatially diversified travel behaviour during non-working days.
Overall, the Local Gi* analysis confirms that spatial clustering of trip activity varies systematically with time, reflecting transitions between commuting and social travel functions. These results demonstrate how temporal segmentation helps reveal distinct behavioural regimes of urban mobility that would otherwise be obscured in aggregated data.
10.5 Mapping \(p\)-Values of Local Gi*
The next analytical step assesses the statistical reliability of the apparent hot and cold zones by mapping the simulated p values corresponding to each \(z\)-score. The mapping function applies fixed classification intervals to highlight confidence levels at 0.001, 0.01, 0.05, 0.10, and 1.00. Lower \(p\) values appear in darker shades, identifying regions where clustering is statistically significant, whereas lighter tones mark areas consistent with random variation. The visual format mirrors that of the z-score maps to enable direct cross-comparison. By reviewing both the intensity and the p value maps together, the analyst can separate visually strong but insignificant clusters from genuinely significant hot or cold areas. This process enhances analytical reliability by ensuring that observed spatial concentrations are statistically defensible rather than products of chance. The mapped \(p\) values therefore serve as the validation layer preceding categorical classification.
The \(p\)-value maps for the Local Gi* analysis provide statistical confirmation of where spatial clustering patterns are most reliable across different time periods. During the weekday morning period, dark red hexagons dominate the western, northeastern, and central sectors of the island, indicating \(p\)-values below 0.01. These areas reflect highly significant clusters of trip activity, suggesting consistent morning travel concentration within residential and employment zones. The dense grouping of low \(p\)-values confirms that the observed spatial patterns are unlikely to occur by random chance and represent genuine hot spot formations.
In the weekday afternoon period, significant areas become more dispersed, and the total number of hexagons with \(p\)-values below 0.01 slightly decreases. This reduction implies greater variability in afternoon trip behaviour as travel flows diversify between work return journeys and secondary activities. Nonetheless, certain stable clusters persist, suggesting enduring transport corridors and mixed-use nodes with high activity density.
The weekend and holiday map shows renewed intensification of significant areas, particularly in central and eastern regions. The spread of dark red hexagons indicates broader spatial clustering of leisure-related movement. This distribution implies that recreational and shopping trips become the dominant source of spatial dependence during non-working days.
Across all three time windows, the persistence of low \(p\)-values validates the robustness of the Local Gi* results. The consistent appearance of significant cells across multiple periods demonstrates that spatial clustering in urban mobility is a stable feature of the transport landscape rather than a transient phenomenon driven by temporal fluctuations.
10.6 Hot/Cold spot classification at \(p\) < 0.05
The final stage consolidates all previous computations and maps into a categorical representation of statistically significant clusters. Using the mapping function map_cluster, each hexagon is classified into one of three categories: Hot Spot, Cold Spot, or Insignificant. The classification depends on both the direction of the Gi* \(z\)-score and the associated p value threshold of less than 0.05. The colour scheme is fixed, with red indicating Hot Spots, blue indicating Cold Spots, and grey indicating Insignificant cells. This uniform convention allows easy visual comparison between different temporal windows. The resulting maps provide a concise yet comprehensive summary of spatial clustering across the study area, clearly identifying where trip activity intensifies or weakens over time. These classified outcomes represent the most interpretable stage of the Gi* workflow, forming the empirical basis for subsequent discussion and interpretation of spatial mobility behaviour.
# --- Classify Gi* into clusters (α = 0.05)--------------------------------------------# p <= .05 significant; Hot if Gi>0, Cold if Gi<0; else InsignificantHCSA_classify <-function(sf_obj, alpha =0.05) { sf_obj |> dplyr::mutate(HCSA_cluster = dplyr::case_when(is.finite(p_sim) & p_sim <= alpha &is.finite(Gi) & Gi >0~"Hot spot",is.finite(p_sim) & p_sim <= alpha &is.finite(Gi) & Gi <0~"Cold spot",is.finite(p_sim) & p_sim > alpha ~"Insignificant",TRUE~NA_character_ ),# lock the factor order to control legend orderHCSA_cluster =factor(HCSA_cluster,levels =c("Insignificant","Hot spot","Cold spot"),ordered =TRUE) )}HCSA_am_c <-HCSA_classify(HCSA_am)HCSA_pm_c <-HCSA_classify(HCSA_pm)HCSA_wk_c <-HCSA_classify(HCSA_wk)# --- Plot helper (fixed order + explicit colors, no “Missing”) -----------------------cluster_levels <-c("Insignificant","Hot spot","Cold spot")cluster_colors <-c("grey80","red","blue") # same order as levelsmap_cluster <-function(sf_obj, title_text) {# drop NA so the legend doesn’t show “Missing” sf_plot <- sf_obj |> dplyr::filter(!is.na(HCSA_cluster)) |># re-assert factor levels in case something changed upstream dplyr::mutate(HCSA_cluster =factor(as.character(HCSA_cluster),levels = cluster_levels,ordered =TRUE))tm_shape(sg_outline) +tm_borders(col ="grey40", lwd =0.7) +tm_shape(sf_plot) +tm_fill(col ="HCSA_cluster",palette = cluster_colors, # categorical palette (tmap v4)showNA =FALSE,title ="HCSA Cluster") +tm_borders(col ="grey85", lwd =0.4) +tm_layout(main.title = title_text,main.title.position="center",main.title.fontface="bold",legend.position =c("right","bottom"),legend.frame =TRUE )}# --- Draw the three maps --------------------------------------------------------------cm_am <-map_cluster(HCSA_am_c, "HCSA Clusters (p < 0.05) — Morning (06–09)")cm_pm <-map_cluster(HCSA_pm_c, "HCSA Clusters (p < 0.05) — Afternoon (17–20)")cm_wk <-map_cluster(HCSA_wk_c, "HCSA Clusters (p < 0.05) — Weekend / Holiday (11–20)")cm_am; cm_pm; cm_wk
The Hot and Cold Spot Analysis (HCSA) cluster maps for weekday morning, weekday afternoon, and weekend or holiday periods collectively reveal the evolution of Singapore’s spatial mobility intensity across temporal contexts, each analysed at \(p\) less than 0.05 significance. In the weekday morning map, a dense network of red hexagonal clusters signifies statistically significant hot spots of trip activity, particularly concentrated across the northern(Yishun, and Woodland), Northern East (Punggol, Sengkang, Hougang), Western (Jurong West, Bukit Batok), Central (Tao Payoh) and Eastern (Tampines, Bedok) portions of the island. These clusters represent high trip-generation zones that likely correspond to large-scale residential districts where outbound commuter traffic peaks between 6 to 9 a.m. The absence of blue cold spot cells indicates that low trip areas during this period follow a more random pattern without statistical clustering. The observed structure suggests a strong and spatially dependent morning mobility framework aligned with commuting towards central employment hubs.
During the weekday afternoon period from 17 to 20 hours, the cluster distribution remains extensive but exhibits slightly lower density and increased dispersion. The presence of multiple smaller red patches indicates a diffusion of activity as people travel home or make short inter-district trips after work. This spatial fragmentation signifies a shift from collective commuter movements to more individualised trip purposes such as school pick-ups, errands, or leisure engagements. The patterns reinforce the dynamic transformation of Singapore’s travel behaviour from highly centralised in the morning to decentralised in the evening.
By contrast, the weekend or holiday map displays broader and more evenly distributed red clusters across central, southern, and eastern zones. These results reflect the transition from occupational to recreational mobility, as travel demand is driven by social, retail, and leisure activities. Unlike weekdays, trip concentration expands into non-core areas, highlighting the more balanced nature of urban activity during off-work periods.
In practical terms, these spatial patterns emphasise how weekday commuting follows structured residential-employment flows, whereas weekend mobility reflects lifestyle-oriented dispersion. Transport planners can use these temporal signatures to synchronise service frequency and route allocation, ensuring capacity aligns with empirically verified high-demand corridors during peak weekday hours and wider area accessibility on weekends.
11 Emerging Hot Spot Analysis (EHSA)
To set the stage, this section introduces a space time approach that asks a practical question about neighbourhood bus mobility. Where are statistically meaningful clusters of higher or lower origin trips appearing, disappearing, or changing through time, and how should we label them in a way that helps planning and policy work. We construct the analysis around the Getis Ord Gi* local statistic that is evaluated for each cell in each time bin and then we examine the sequence of those statistics using the Mann Kendall trend test. The combined outcome is a categorical label such as new hot spot, intensifying hot spot, persistent hot spot, diminishing hot spot, or their cold analogues. We run the procedure for three mandated periods that reflect commuting and leisure rhythms, and we carry the results back to the validated mainland hexagon grid so that interpretation is always spatially anchored. The purpose is not simply to produce colourful maps but to generate defensible evidence about evolving trip concentration and suppression across the urban fabric.
11.1 Build neighbours and weights on hex grid
At the core of any local cluster method is a clear statement of spatial relationships. This subsection prepares that foundation by defining which hexagons are considered neighbours and how much influence each neighbour should exert in the calculations. We adopt contiguity based neighbours with the queen rule so touching at an edge or a corner counts, and we include each cell itself to ensure that isolated or island cells can still be processed without dropping from the analysis. We then compute inverse distance weights from the centroids so that very close neighbours have more influence than distant ones while still keeping the structure simple and reproducible. Warnings about cells with no neighbours can appear in coastal or fragmented areas, and the include self step addresses that gracefully. Preparing neighbours and weights on the geometry first and storing them as named columns allows the later space time cube to attach the correct structure without costly recomputation at every time step.
# ---- 11.0 Build neighbours & weights on hex grid -----------------------# Assumes we already have:# - trips_panel_sf: columns HEX_ID, DAY_TYPE, HOUR_OF_DAY, TRIPS (sf ok; we will drop geometry)# - hexagon_active: sf polygons with HEX_ID + geometry# - sg_outline: sf polygon of SG mainland (optional, for basemap)hexagon_geo <- hexagon_active |>mutate(nb =include_self(st_contiguity(geometry, queen =TRUE)),wt =st_inverse_distance(nb, geometry, scale =1, alpha =1) )# NOTE: "some observations have no neighbours" warnings are expected for islands;# include_self() ensures each hex has at least itself.
11.2 Period definitions
In this step we translate the assignment brief into explicit time windows that a computer will understand. The weekday morning period covers the early departure window from 06:00 to 09:00 that captures school and work origins, the weekday afternoon period covers the return window from 17:00 to 20:00 that mixes home bound and after work activities, and the weekend or holiday period focuses on late morning through early evening from 11:00 to 20:00 when leisure and retail activity dominate. Declaring these ranges once, with a readable label for each, prevents accidental drift in later filters and guarantees that any figures or tables can be regenerated consistently. The labels also become facets in charts and map titles, which helps a reader keep track of the storyline. Choosing hour as the time bin aligns with the data structure of origin trips and keeps the sample size per bin adequate for reliable local statistics and trend testing.
Before running the main classification, we safeguard the workflow against differences in package versions. The emerging hotspot procedure returns the trend result either as a simple column named sl, as a column named mk_sl, or nested inside a list column that contains the Mann Kendall (MK) object. The helper in this subsection extracts the \(p\) value reliably no matter which structure appears, and falls back to other common names if needed. Returning a numeric vector rather than a mixed object means that downstream filtering is straightforward and readable. This design choice protects the analysis from breaking on another machine or in a future semester when packages have changed. It also makes the significance threshold explicit and therefore auditable. By separating this concern into a short utility we keep the runner function focused on analysis logic, while confidence in reproducibility and portability is improved for both code review and grading.
# ---- helper to pull MK p-value robustly --------------------------------.get_mk_p <-function(tbl) {# Return a numeric vector pval; works for multiple sfdep versions nms <-names(tbl)if ("sl"%in% nms) return(tbl$sl)if ("mk_sl"%in% nms) return(tbl$mk_sl)if ("mk"%in% nms) {return(purrr::map_dbl(tbl$mk, function(x) {if (!is.null(x$sl)) return(as.numeric(x$sl))if (!is.null(x$p.value)) return(as.numeric(x$p.value))if (!is.null(x$pvalue)) return(as.numeric(x$pvalue))NA_real_ })) }# Fallback: look for generic p-value fields if presentif ("p_value"%in% nms) return(tbl$p_value)if ("pvalue"%in% nms) return(tbl$pvalue)if ("p.val"%in% nms) return(tbl$`p.val`)rep(NA_real_, nrow(tbl))}
11.4 EHSA runner
In this section we operationalise the full analytical sequence for a single period. The function accepts a day type, a set of hour values, and a human friendly label. It filters the balanced panel to exactly those rows, constructs a space time cube keyed by hexagon identifier and hour, and attaches the previously prepared neighbours and weights using their column names. With this structure in place, the emerging hotspot procedure computes a series of Gi star statistics per cell across the chosen hours, performs the Mann Kendall test on that temporal sequence, and assigns a classification label that captures both significance and direction of change. We then keep only features with significant trend evidence, and we join the result back to the hexagon geometry so that mapping and counting become trivial. Encapsulating all of this in one function makes the workflow concise, testable, and repeatable across periods without copy paste mistakes.
Next we scale the runner across the three mandated periods in a single call using a functional pattern. A parameter map sends the weekday morning, weekday afternoon, and weekend or holiday arguments into the runner, gathers the three outputs, and merges them into one tidy table with a period column. This structure supports side by side cartography, compact summary tables, and faceted charts without additional reshaping. The approach reduces the risk of silent inconsistencies that can arise when running separate code blocks and also speeds iteration when we refine weights or significance thresholds. If new periods are added in future work, they can be introduced by editing only the small period table rather than rewriting logic. The result is a clean separation of configuration and execution that supports both transparency for assessment and flexibility for further inquiry.
# ---- Run for all periods ------------------------------------------ehsa_list <- purrr::pmap(list(periods$day_type, periods$hours, periods$label), run_ehsa)ehsa_all <- dplyr::bind_rows(ehsa_list)
11.6 Quick class distribution
Before drawing maps, it is useful to perform a numerical sense check that confirms the outcome is plausible. The bar chart in this section counts the number of hexagons in each classification for every period and presents them in small multiples. This quick view reveals whether the analysis is returning only a few scattered features or a broad swath of classifications, and it highlights imbalances such as many cold labels and very few hot labels or the opposite. Extreme skews can suggest parameter issues, for example an overly sparse neighbour structure or an hour range that is too narrow. The figure also provides the quantitative backbone for the narrative paragraphs in the report, because we can cite counts while the maps supply the geography. Keeping this diagnostic small and fast encourages repeated use whenever we adjust upstream choices.
# ---- 11.4 Quick class distribution ---------------------------ggplot(ehsa_all, aes(x = classification)) +geom_bar() +facet_wrap(~ period, scales ="free_y") +labs(title ="EHSA classes (significant only, p < 0.05)",x =NULL, y ="Hexagon count") +theme_bw() +theme(plot.title =element_text(hjust =0.5, face ="bold"))
The bar chart tells a clear story about temporal regimes. In the weekday morning panel and the weekday afternoon panel, almost every significant cell falls into a single class while the counts for all other classes are negligible. In practice this reads as a broad field of persistent cold conditions during the short commuting windows. Local Gi star values at those hours do not trend upward through time and the Mann Kendall test confirms that the sequences are consistently low rather than moving toward higher intensity. This aligns with a mature and stable commuting pattern where origins are widely distributed and rarely form accelerating clusters within a three or four hour slice. The immediate takeaway for operations is to protect reliability rather than pursue large reallocations. Priority at known choke points, disciplined headway control, and dwell time management will matter more than changing routing inside these windows. Monitoring is still necessary, but the evidence does not call for structural changes on weekdays.
The weekend or holiday panel looks very different. Counts spread across several classes and no single label monopolises the distribution. The largest bar is a cold class that reflects enduring low origin intensity across green and water dominated areas and in some employment zones that quieten on non working days. At the same time there are real pockets of demand growth. Intensifying hot classes and consecutive hot classes indicate cells where sequences of local Gi star statistics have moved upward across the late morning to evening bins. These cells typically coincide with town centres, interchange nodes, waterfront leisure corridors, and retail districts. Sporadic hot classes appear in places that likely host irregular events, suggesting transient surges that are worth watching but may not justify permanent capacity.
Planning implications are direct. Weekday peaks show steadiness, so keep the focus on reliability and enforcement of service levels, with small tactical boosts only where crowding is observed. Weekend daytime deserves targeted investment. Add trips or deploy larger vehicles on corridors that connect residential belts to regional centres and recreation zones. Use the high count classes to prioritise stops for queue management and passenger information. Finally, treat these results as a screening layer. Re run the analysis with a knn neighbour scheme as a sensitivity check, extend the weekend window to confirm persistence, and compare outcomes with observed load factors to translate spatial statistics into service adjustments that riders will feel.
Note
We can use below code chunk to generate a compact frequency table of EHSA classes by period. It groups by period and classification, counts hexagons in each combination, and sorts from most common to least. The table supports quick plausibility checks, helps identify dominant hot or cold patterns, and provides precise numbers to reference in the narrative and figure captions.
The frequency table confirms what the maps suggested. During the weekday morning period and the weekday afternoon period, almost every hexagon carries no significant emerging signal. Each panel reports about 826 cells with missing class values, which in this workflow reflects features that failed the significance screen and therefore received no label after the join back to geometry. In other words, within the short commuting windows the local Gi star sequences are largely stable rather than trending upward or downward, so the Mann Kendall test does not indicate meaningful change. This supports an operations focus on reliability and crowd management rather than network redesign inside those hours.
The weekend or holiday period shows greater variety, although a large share still carries no significant change with about 433 cells. Among the labelled results, the most common class is no pattern detected with about 245 cells, followed by sporadic hotspot with about 41 cells and persistent coldspot with about 39 cells. Sporadic coldspot appears about 36 times, intensifying hotspot about 19 times, consecutive hotspot about 7 times, and new coldspot appears a few times. The mix tells us that late morning to evening on non working days has real but spatially contained growth around town centres, leisure corridors, and interchange areas, while many other places remain consistently quiet. Service planning should therefore prioritise targeted weekend boosts on the corridors and stops that fall in intensifying or consecutive hotspots, while weekday peaks should emphasise headway control, signal priority, and dwell discipline to keep established flows moving smoothly.
11.7 Visualising EHSA maps in different windows
To communicate results clearly to readers who first look at plots, this section defines a single mapping function that produces publication quality thematic maps with consistent design choices. The function adds a neutral mainland outline for context, fills hexagons by the classification field using a readable categorical palette, draws light borders to preserve cell shapes, and places the legend outside the frame so that the data area stays uncluttered.
plot_ehsa_map(filter(ehsa_all, period =="Weekday Afternoon Peak"),"EHSA of Origin Trips — Weekday Afternoon Peak")
plot_ehsa_map(filter(ehsa_all, period =="Weekend/Holiday Peak"),"EHSA of Origin Trips — Weekend or Holiday Peak")
Reading the three EHSA maps together shows a sharp contrast between weekday commuting periods and weekend activity. In the weekday morning map almost all analytical cells are white, which means the Gi star sequences across the hours from six to nine do not exhibit a consistent upward or downward tendency that passes the Mann Kendall test at the five percent level. This does not mean there is no morning demand. Rather it signals that within this tight three hour window the local intensity is stable or noisy instead of directional. Where a few isolated coloured cells appear they do not form a coherent structure across the island and therefore carry little strategic meaning on their own.
A very similar message appears in the weekday afternoon map. Again most cells are not significant. The period from seventeen to twenty captures home bound and after work travel that is well established in the network. The absence of emerging hot or cold classes suggests a mature pattern with reliable peaks but without strong within period trajectories. For operations this favours actions that protect reliability, such as priority measures, dwell control, and headway management, more than changes in routing or capacity distribution inside the short window.
The weekend or holiday map tells a different story. Many coloured cells appear and several categories are visible at once. Persistent cold classes line up with water bodies and large parks where bus origins remain low on weekends, while new or sporadic cold classes suggest employment districts that quieten when offices close. In contrast, persistent and intensifying hot classes cluster around major town centres, interchange nodes, and waterfront leisure corridors where shopping and recreation attract sustained flows from late morning to evening. Consecutive hot classes near dense transfer areas point to a recent pickup that persisted into the last hour set, and sporadic hot classes around event focused precincts hint at irregular surges.
The practical implications are clear. Weekday peaks look structurally stable, so focus resources on keeping services reliable and clearing bottlenecks. Weekend midday to evening shows genuine growth corridors, so plan extra capacity on routes that connect residential belts to regional centres, parks, and coastal attractions. Monitor cells that moved to intensifying or consecutive hot classes because they are strong candidates for targeted frequency increases, stop level crowd management, or demand responsive extras on specific corridors. Finally, remember that the trend test uses short hour series, so it is worth validating these findings with a longer weekend window or daily bins and by checking an alternative neighbour structure to confirm robustness.
12 Discussion
Here’s a direct, evidence-based answer to the 4 research questions, grounded in this study and its EHSA/LMSA workflow.
RQ1 — Spatial distribution
Weekday demand is strongly diurnal and spatially patterned. Our own peak-hour maps show that weekday mornings (06:00–09:00) concentrate origins in northern and eastern heartlands—Woodlands, Yishun, Sengkang, Tampines, Bedok, Bukit Panjang and Jurong West—consistent with home-to-work flows leaving residential areas for employment centres. By contrast, weekday afternoons (17:00–20:00) shift intensity toward the Central Region (Orchard, Downtown Core, Kallang), reflecting return trips and after-work activity. Weekends and holidays (11:00–20:00) exhibit a broader, flatter profile with dispersed origins that remain prominent around retail and leisure corridors and town centres from late morning into evening. These patterns align with the diagnostic hourly profile that shows two clear weekday peaks and a single broad weekend peak, validating the period definitions used in the analysis.
RQ2 — Local spatial clustering
The study applies Local Moran’s I, and Getis-Ord Gi* to detect where high or low ridership clusters and spatial outliers occur, moving beyond global summaries to location-specific evidence. These local statistics reveal neighbourhoods of consistently high or low values and pinpoint outliers that diverge from their surroundings, strengthening the descriptive reading of the peak maps. In practice, the observed morning heartland concentrations and afternoon central-area intensifications manifest as high-value clusters (hot clusters) in residential belts during AM and around core employment/retail zones during PM, while water-dominated or park areas and certain weekend-quiet employment precincts register as low-value clusters (cold clusters). Using LMSA in tandem with the period-specific mapping thus ties intuitive patterns to statistically significant local structure, ensuring that subsequent policy discussion rests on formal evidence rather than visual impression alone.
RQ3 — Temporal evolution (EHSA with Mann–Kendall)
EHSA on space–time Local Getis-Ord Gi* sequences shows a clear contrast between weekday peaks and weekends. For weekday morning and afternoon windows, most hexagons are not significant in the Mann–Kendall trend test; the short three-hour windows appear stable or noisy rather than trending upward or downward, indicating mature, well-established commuting peaks without strong within-period trajectories. On weekends/holidays, many cells become significant and a variety of EHSA classes emerge. Persistent and intensifying hot spots cluster around town centres, interchange nodes and waterfront leisure corridors, signalling sustained or strengthening midday-to-evening flows. Persistent cold spots coincide with water bodies and large parks; “new” or sporadic cold spots appear in office districts that quieten when workplaces are closed. Consecutive and sporadic hot spots around event-oriented precincts indicate recent or irregular surges that merit monitoring. Together, these results show stability on workday peaks but genuine growth and variability on weekends.
RQ4 — Policy and planning implications
Because weekday peaks are structurally stable within their short windows, operations should prioritise reliability over major structural changes: protect headways, manage dwell times, and address recurrent bottlenecks rather than reallocating capacity inside the three-hour peaks. In contrast, weekend/holiday midday–evening corridors with persistent or intensifying hot spots warrant targeted capacity increases—additional trips or larger vehicles—especially on links from residential belts to regional centres, parks and coastal attractions. Cells that moved into intensifying or consecutive hot classes are prime candidates for stop-level crowd management, queue design, and real-time passenger information. As a good practice, treat EHSA as a screening layer: validate hot/cold trends against observed loads, test neighbour definitions (e.g., k-nearest-neighbours sensitivity), and consider extending the weekend window to confirm persistence before committing long-term resources. These steps translate the spatial–temporal evidence into actionable network adjustments that riders will feel.
13 Conclusion
This study demonstrates that a careful combination of credible data, defensible spatial units, and rigorous local statistics can turn raw ridership counts into actionable intelligence for a bus network. Using LTA DataMall origins linked to the mainland bus-stop system and aggregated to a validated 750-metre hexagon grid, we built a fully balanced space–time panel and applied two complementary lenses: local measures of spatial autocorrelation to locate meaningful clusters, and EHSA to track how those clusters evolve within defined peak windows. The workflow—projection to SVY21, mainland masking, active-cell filtering, stable IDs, harmonised cartography, neighbour and weight specification, and reproducible code—ensures that every figure and metric can be regenerated and audited. 
Three broad findings follow. First, weekday peaks are structurally stable. Both morning and afternoon windows show very few cells with significant emerging trends, indicating mature commuting patterns that repeat predictably within those tight three-hour slices. This does not imply weak demand; rather, it shows that within-window Gi* sequences are steady enough that the Mann–Kendall test rarely flags change. Second, weekends and holidays are different. EHSA returns a diverse mix of classes, with persistent and intensifying hot spots around town centres, major interchanges, and waterfront leisure corridors, while persistent cold spots coincide with parks, reservoirs, and water. Sporadic and consecutive hot spots near event-oriented districts reveal irregular but important surges that deserve operational attention. Third, the spatial story seen in descriptive maps is confirmed statistically: morning origins concentrate in northern and eastern heartlands, afternoon activity tilts toward the Central Region, and weekend flows broaden along retail and recreation axes.
Policy implications are direct. On weekdays, the priority should be reliability: protect headways, manage dwell times, relieve recurring choke points, and use targeted bus-priority or stop-level design to keep predictable peaks moving. On weekends and holidays, deploy flexible capacity where the evidence shows growth—extra trips or larger vehicles on links between residential belts and leisure or retail clusters; queue and crowd management at stops inside intensifying or consecutive hot spots; and real-time information to smooth surges. Finally, treat EHSA as a screening tool and institutionalise sensitivity checks: test alternative neighbour schemes, extend or shift weekend windows when warranted, and triangulate with observed loads. Taken together, these steps translate spatial analytics into service choices that improve rider experience while using resources where they matter most.
14 References
Boschan, J.A. and Roman, C.G. (2024) ‘Hot spots of gun violence in the era of focused deterrence: A space-time analysis of shootings in South Philadelphia’, Social Sciences, Vol. 13.
Kam, T.S. (2025) ‘Chapter 8: Spatial Weights and Applications’, R for Geospatial Data Science and Analytics (R4GDSA). Available at: https://r4gdsa.netlify.app/chap08.html (Accessed: October 2025).
Kam, T.S. (2025) ‘Chapter 9: Global Measures of Spatial Autocorrelation’, R for Geospatial Data Science and Analytics (R4GDSA). Available at: https://r4gdsa.netlify.app/chap09.html (Accessed: October 2025).
Kam, T.S. (2025) ‘Chapter 10: Local Measures of Spatial Autocorrelation’, R for Geospatial Data Science and Analytics (R4GDSA). Available at: https://r4gdsa.netlify.app/chap10.html (Accessed: 23 October 2025).
Kam, T.S. (2025) ‘In-Class Exercise 06: Emerging Hot Spot Analysis (EHSA)’, ISSS626 Geospatial Analytics and Applications (AY2025/26 – August Term). Available at: https://isss626-ay2025-26aug.netlify.app/in-class_ex/in-class_ex06/in-class_ex06-ehsa#/title-slide (Accessed: October 2025).
Kim, M. and Lee, S. (2023) ‘Identification of emerging roadkill hotspots on Korean expressways using space–time cubes’, International Journal of Environmental Research and Public Health, 20(6).
Kam, T.S. (2025) ‘Lesson 4: Spatial Weights and Applications’, ISSS626 Geospatial Analytics and Applications (AY2025/26 – August Term). Available at: https://isss626-ay2025-26aug.netlify.app/lesson/lesson04/lesson04-spatial_weights#/title-slide (Accessed: October 2025).
Kam, T.S. (2025) ‘Lesson 5: Global Spatial Autocorrelation (GLSA)’, ISSS626 Geospatial Analytics and Applications (AY2025/26 – August Term). Available at: https://isss626-ay2025-26aug.netlify.app/lesson/lesson05/lesson05-glsa#/title-slide (Accessed: October 2025).
Mack, Z.W.V. and Kam, T.S. (2018) ‘Is there space for violence?: A data-driven approach to the exploration of spatial-temporal dimensions of conflict’, Proceedings of the 2nd ACM SIGSPATIAL Workshop on Geospatial Humanities, Seattle, WA, USA, 6 November 2018, pp. 1–10. ACM Digital Library. doi: 10.1145/3282933.3282935.
Tan, Y.Y. and Kam, T.S. (2019) ‘Exploring and visualizing household electricity consumption patterns in Singapore: A geospatial analytics approach’, Information in Contemporary Society: 14th International Conference, iConference 2019, Washington, DC, March 31–April 3, 2019: Proceedings, pp. 785–796. Springer. doi: 10.1007/978-3-030-15742-5_74.
Transportation Modelling Group, University of Toronto (2021) Traffic Zone Guidance: March 2021 (Final Report). Toronto: University of Toronto. Available at: https://tmg.utoronto.ca/files/Reports/Traffic-Zone-Guidance_March-2021_Final.pdf (Accessed: October 2025).
